We start with the definition of the problem: given a set 'l' and a
number 'n', return the set of all subsets of 'l' of cardinality
'n'. Sets are represented by lists.

The key was to choose the right definition.

Let ps(L) is a powerset of L that does not include {empty set} (that
is, a singular set whose element is the empty set).
        length(ps(L)) = (- (expt 2 (length L)) 1)

Let L be a non-empty set and let L = A U B where A and B are two
disjoint subsets of L.

Obviously,
        ps(L) = ps(A U B) =
        ps(A) U ps(B) U { y U x | x <- ps(A), y <- ps(B) }
Let (subsets L n) = (filter (lambda (el) (= n (length el))) ps(L) )

Thus the desired function 'subsets' is a filtered powerset ps. This
seems to be a stupid definition. Let's not be hasty however. Note that
filter and union commute: the filter of a union of sets is the union
of filtered sets. Therefore, from the previous expression

        (subsets L n) = (subsets (union A B) n)
        = (subsets A n) U (subsets B n)
        Union{ y U x | x <- (subsets A k), y <- (subsets B (- n k)),
               k=1,n-1 }

Well, this is it. Note, we didn't say how to split L into two disjoint
subsets A and B. We can do as we wish. For example, we can choose to
split in such a way so that (length B) is n. In the most difficult
case where n = (/ (length L) 2), this corresponds to a "divide and
conquer" strategy, so to speak (at least in the first stages).

The Scheme code below implements this idea verbatim. It uses the
accumulator-passing style that was expounded earlier.

(define (subsets-v5 l n)

  ; The initialization function. Check the boundary conditions
  (define (loop l ln n accum)
    (cond
     ((<= n 0) (cons '() accum))
     ((< ln n) accum)
     ((= ln n) (cons l accum))
     ((= n 1)
      (let fold ((l l) (accum accum))
        (if (null? l) accum
            (fold (cdr l) (cons (cons (car l) '()) accum)))))
     (else
      (split l ln n accum))))

  ; split l in two parts a and b so that (length b) is n
  ; Invariant: (equal? (append a b) l)
  ; ln is the length of l
  (define (split l ln n accum)
    (let loop  ((a '()) (b l) (lna 0) (lnb ln))
      (if (= lnb n) (cont a lna b lnb n accum)
          (loop (cons (car b) a) (cdr b) (+ 1 lna) (- lnb 1)))))

  ; The main body of the algorithm
  (define (cont a lna b lnb n accum)
    (let* ((accum
            (loop a lna n accum))
           (accum       ; this is actually (loop b lnb n accum)
            (cons b accum))
           )
      (let inner ((k 1) (accum accum))
        (if (> k (min lna (- n 1))) ; don't loop past meaningful boundaries
            accum
            (let ((as (loop a lna k '()))
                  (bs (loop b lnb (- n k) '())))
              (inner (+ 1 k)
                     ; compute the cross-product of as and bs
                     ; and append it to accum
                     (let fold ((bs bs) (accum accum))
                       (if (null? bs) accum
                           (fold (cdr bs)
                                 (append
                                  (map (lambda (lst) (append lst (car bs))) as)
                                  accum))))))))))

  (loop l (length l) n '()))

The benchmark runs in 993 ms of user time and allocates only 36.5 MB
of memory, on Gambit-C interpreter. This is the absolute, incredible
record. Under SCM:

subsets-v3 (called combos by John David Stone)
;Evaluation took 1596 mSec (98 in gc) 657662 cells work, 4721364 env, 97 bytes other

subsets-v5:
;Evaluation took 700 mSec (322 in gc) 1708112 cells work, 610264 env, 105 bytes
other                                                                          
That is, more than twice as fast.

Continuing the table from the previous post:

Procedure       Gambit-C               Bigloo 2.4b      Bigloo 2.4b
                interpreter, s         interpreter, s   compiler, s
subsets-v0        285.0                   11.59          5.62    3.14
subsets-v1          6.3                    5.45          2.22    0.34
subsets-v3          8.1                    4.78          0.96    0.27
subsets-v20        14.1                    5.53          0.96    0.26    
subsets-v21         7.7                    4.88          0.66    0.26
subsets-v22         5.0                    3.18          0.62    0.25
subsets-v23         4.1                    2.86          0.82    0.25
subsets-v5          0.9                    1.56          1.10    0.76

Well, compiled code isn't that fast -- but that because subsets-v5
isn't too optimal. The append and map ought to be converted into the
accumulation-passing style. That will reduce the amount of garbage as
well. But it's past midnight.

    But what _is_ subsets function?  What does it
compute, and what is it used for?

    I would be grateful for additional comments.

        That is:  Apart from garbage collection, the two procedures are
very closely comparable -- but SUBSETS-V3 creates only a third as much
garbage.

        My guess is that, in my testing environment, SUBSETS-V3 is getting
a bigger advantage from Chez Scheme's optimizations.  The results I get under
MzScheme are similar to Oleg's:

And he's right.

  http://sardine.cs.columbia.edu:8080/subsets/

All graphs are timings of different lengths out of a set of 22
elements, msubset is the memoization of Oleg's first version, msubset1
is even simpler, removing the special case for len=1.  The code is at

  http://sardine.cs.columbia.edu:8080/subsets/subsets.ss

but this is just for reference, it contains all the code pieces with
some junk that makes gnuplot doo all the graphs.

Obviously a lazy functional language like Haskell shows this
implementation in its best light.  As Oleg posted, we can define the
subsets function mathematically as

subsets(S, 0)       = {{}}
subsets({}, n)      = {}
subsets({a} U S, n) = {{a}} x subsets(S, n-1)  U  subsets(S, n)

This translates directly into Scheme or Haskell.  Scheme code has
already been posted.  In Haskell it looks like this:

subsets _      0 = [[]]
subsets []     _ = []
subsets (x:xs) n = map (x :) (subsets xs (n - 1)) ++ (subsets xs n)

This is slow, partly because it conses too much but mostly because the
same subproblems are solved repeatedly.  Memoization can be used to
attack this performance problem.  The lazy streams approach will give
many the benefits of memoization while avoiding most of the
disadvantages.

Memoization pros

 - can be applied automatically without changing the naive
   implementation

 - memoized results can be saved across top level invocations of the
   function

Memoization cons

 - not always faster

 - doesn't work well if function arguments are complex (the required
   hash compare might be too slow)

 - usually implemented using imperative machinery under the hood
   (although the resulting memoized function has functional behavior)

 - can lead to horrendous space leaks in the memo tables; recovering
   this space might require manual control partly negating the first
   pro listed above

In a lazy stream formulation, the nth element of the lazy stream will
be a list of the subsets of cardinality n.  In other words, the stream
will be [subsets(S,0), subsets(S,1), subsets(S,2), ...].  The first
stream element will always be the list of the empty set.  For all n
greater than length of the input list the element will be null.  This
is very smooth in Haskell:

subsets s n = (subsets_stream s) !! n

subsets_stream []     = [[]] : repeat []
subsets_stream (x:xs) =
    let r = subsets_stream xs
        s = map (map (x:)) r
    in [[]] : zipWith (++) s (tail r)

As a small bonus We can easily use subsets_stream to implement the
powerset function:

powerset s = concat (takeWhile (not . null) (subsets s))

Unfortunately this isn't as pretty in Scheme.  Comparing the Scheme
implementation below to Haskell, Scheme is far more verbose since we
have to explicitly request laziness, we have to provide lazy
equivalents of standard library functions (CONS, CAR, CDR, LIST-REF
and MAP, although these could be put in an external library) and
finally currying is more convenient to use in Haskell than it is in
Scheme.  (Often Haskell's pattern matching would be an advantage as
well, but here it plays no important part.)  One tricky part about
explicitly requesting laziness is that it's easy to get wrong.  Look
at the mix of uses of STREAM-MAP vs MAP and STREAM-CONS vs CONS in
SUBSET-STREAMS below.

Since most Scheme implementations don't natively support R5RS macros,
I have simply replaced all uses of STREAM-CONS with
(CONS (DELAY a) (DELAY b)).  Notice that I use fully lazy streams.
The streams in SICP have lazy tails but eager heads.

; Most Scheme implementations don't have R5RS macros, alas.
;(define-syntax stream-cons
;  (syntax-rules ()
;    ((stream-cons hd tl) (cons (delay hd) (delay tl)))))

(define (subsets s n)
  (stream-ref (subsets-stream s) n))

;; Return a lazy stream of all subsets of list XS.  The nth element of
;; the stream is a list of the subsets of cardinality n (counting from 0).
(define (subsets-stream xs)
  (if (null? xs)
      (cons (delay (list '())) (delay (stream-repeat '())))
      (let* ((r (subsets-stream (cdr xs)))
             (s (stream-map (lambda (ss)
                              (map (lambda (y) (cons (car xs) y)) ss))
                            r)))
        (cons (delay (list '()))
              (delay (stream-map append s (stream-cdr r)))))))

;; The remainder should be part of a lazy streams library.
;; These are the lazy list equivalents of CAR, CDR, LIST-REF and MAP.
;; STREAM-REPEAT is inspired by Haskell's repeat function.

(define (stream-car s)
  (force (car s)))

(define (stream-cdr s)
  (force (cdr s)))

;; Create an infinite stream of X's.
(define (stream-repeat x)
  (cons (delay x) (delay (stream-repeat x))))

(define (stream-ref s n)
  (if (= n 0)
      (stream-car s)
      (stream-ref (stream-cdr s) (- n 1))))

;; Answer to an easy exercise in SICP.  Most SICP exercises are harder ;-)
(define (stream-map f . streams)
  (if (null? (car streams))
      '()
      (cons (delay (apply f (map stream-car streams)))
            (delay (apply stream-map
                          (cons f (map stream-cdr streams)))))))

  http://sardine.cs.columbia.edu:8080/subsets/

Saving results across top level invocations is at odds with the space
leak problem I mention later.  For a restricted set of arguments for
this specific problem, you can encapsulate the answers for all of the
integer subsets in one place.  This is a bit ugly in Scheme, and it
also changes the arguments to the subsets function.  (Instead of
passing a list and a number we need two numbers.)

(define (ints-from n) (cons (delay n) (delay (ints-from (+ n 1)))))

(define (int-subsets-stream-from n)
  (cons (delay (subset-stream (ints-from 1))
        (delay (int-subsets-stream-from (+ n 1))))))

;; Return a list of all subsets of the integers {1, ..., N} of
;; cardinality R.
(define subsets
  (let ((all-subsets (cons (delay '(())
                           (delay (int-subsets-stream-from 1))))))
    (lambda (n r)
      (stream-ref (stream-ref all-subsets n) r))))

Space leaks are a problem here and would require manual intervention.
(No better than memoization on this point.)

Ignoring the space issues, this is a one-liner in Haskell:

int_subsets = [subsets_stream [1 .. n] | n <- [0 ..]]
subsets n r = (int_subsets !! n) !! r

I didn't mention the performance benefits of memoization from avoiding
repeated computation of the same subproblems since it is an obvious
feature.  I was thinking of it when I wrote "many", but perhaps many
wasn't the right word when I meant just two out of three.

Of course laziness can certainly hurt performance as well.  In Haskell
this problem shows up in the other direction -- you have to decide
when you want to be strict instead of lazy.  My guess is that this is
normally harder than deciding when memoization will help.