We compute the intersection of two subspaces where a subspace is expressed as a list of spanning vectors. Such a list is viewed as a matrix. A dual space is much like an orthogonal space but does not require an inner product. If X is a vector space then X* is the dual of X.

We use the notation that when x is a set of vectors in a vector space then S(x) is the set of vectors that can be expressed as linear combinations of vectors in x. S(x) is the smallest subspace containing x. It is also the intersection of all the subspaces containing x. When p belongs to X*, p is a linear operator mapping X to the reals. The notion of orthogonal space is supplanted by a sub-space in the dual. If Z is a linear subspace of X then Z⊥ = {q∊X*|∀x(x∊X → qx = 0)} is a linear subspace of X*. (X⊥)⊥ = X and for subspaces U and V of X we have (U∩V) =(S(U⊥ ∪ V⊥))⊥.

If {bi} is a basis for X then there is a unique basis {ci} for X* such that cjbi = δij.

If we have a subspace U of X spanned by a set ui of independent vectors, then we will want to find a set of independent vectors in X* which span U⊥; these will not be unique. We will then exploit (U∩V) =(S(U⊥ ∪ V⊥))⊥ to compute the intersection of subspaces from their spanning sets.

Computing a spanning set for U⊥

Here we speak of a k×n matrix spanning a subspace of an n dimensional vector space when the k vectors, each taken from a row of the matrix, together span the subspace. The description “reduced row echelon form” (RREF) refers to a property of a matrix which is loose enough so that any subspace may be spanned some RREF matrix, but also strict enough that just one RREF matrix, except perhaps for some bottom zero rows, spans the subspace. Adding a zero row to a matrix does not change what the matrix spans. For our purposes here we need the dual of RREF. We will say that a matrix is in dual RREF (DRREF) when the result of reversing the order of both the rows and columns of the matrix results in a matrix in RREF.

Suppose we put the vectors that define a subspace U (by spanning) into a matrix by rows and then transform the matrix into reduced row echelon form and call the result R. U determines R uniquely. The row vectors of R also span U. If our original vectors were independent there will be no zero rows in R. We henceforth presume our original vectors to be independent. Some of the columns of R will be all 0 except for one 1 and that 1 will be the first (left most) non zero of its row. We call these columns “sentinels”. Given a matrix R in this form of row vectors that span U, there is a simple process to produce a matrix R* of row vectors that span the U⊥ which is a subspace of X*. The basis for X* will be the dual basis as introduced above. Here are some symmetric relations between R and R* that uniquely determine either from the other:

R and R* each determine the other; the relationship is dual. Here are some examples of R and R* pairs—one above the other:
( 1  3)

(-3  1)
------
(0  1  5)

(1  0  0)
(0 -5  1)
------
( 1  2  0  5)
( 0  0  1  6)

(-2  1  0  0)
(-5  0 -6  1)
------
(0  1  2  6  0  1  0  0  4  0)
(0  0  0  0  1  1  0  0  1  1)
(0  0  0  0  0  0  1  0  4  1)
(0  0  0  0  0  0  0  1  2  1)

(1  0  0  0  0  0  0  0  0  0)
(0 -2  1  0  0  0  0  0  0  0)
(0 -6  0  1  0  0  0  0  0  0)
(0 -1  0  0 -1  1  0  0  0  0)
(0 -4  0  0 -1  0 -4 -2  1  0)
(0  0  0  0 -1  0 -1 -1  0  1)
Perhaps it is clear that the if p ∊ R and q ∊ R* then qp = 0 and so consequently for the spaces spanned for the rows of R and R*.

Here is some Scheme code to do these unnatural matrix acts; (ossp R) returns R*:

; (iden n) yields the n by n identity matrix.
(define (iden k) (if (zero? k) '() (if (= k 1) '((1)) (let ((z (iden (- k 1))))
   (cons (cons 1 (cons 0 (cdar z))) (map (lambda (v) (cons 0 v)) z))))))

; flip negates each element of a matrix.
(define (flip x) (if (null? x) '() (cons 
   (let ifl ((a (car x))) (if (null? a) '() (cons (- (car a)) (ifl (cdr a)))))
   (flip (cdr x)))))

; ((cull index-list) list) deletes numbered elements from list.
; ((cull '(2 4)) '(10 11 12 13 14 15)) => (10 11 13 15)
(define ((cull n) l) (let px ((l l)(n n)(k 0)) (if (null? n) l (if (null? l) '()
    (if (= (car n) k) (px (cdr l) (cdr n) (+ k 1))
      (cons (car l) (px (cdr l) n (+ k 1))))))))

; clz counts leading zeros.
(define (clz a) (if (or (null? a) (not (zero? (car a)))) 0 (+ 1 (clz (cdr a)))))

(define (transpose x) (if (null? (car x)) '()
   (cons (map car x) (transpose (map cdr x)))))

(define (ossp k) (transpose (let* ((a (map clz k)) (b (map (cull a) k)))
  (let puff ((a a)(k (flip b))(i (iden (length (car b))))(j 0))
    (if (and (null? k) (null? i)) '()
       (if (and (not (null? a)) (= j (car a)))
          (cons (car k) (puff (cdr a) (cdr k)      i  (+ j 1)))
          (cons (car i) (puff      a       k  (cdr i) (+ j 1)))
     ))))))
Returning to a previous goal of exploiting (U∩V) =(S(U⊥ ∪ V⊥))⊥ we try the following:
(define (inter a b) (let* ((x (ossp (rref (minspan a))))(y (ossp (rref (minspan b))))
    (z (minspan (append x y)))) (ossp (rref z))))
This works in a few test cases but instead of using the somewhat expensive minspan function it would be better to merely strip the trailing zero vectors from the rref yield. Even better, cause ossp to ignore trailing zero vectors in its argument. With this version oss of ossp:
(define (oss k) (transpose (let* ((a (map clz k))
     (k (let ((l (length (car k))))
     (let ly ((k k) (a a)) (if (null? k) '()
       (if (= l (car a)) '() (cons (car k) (ly (cdr k)(cdr a))))))))
   (b (map (cull a) k)))
   (let puff ((a a)(k (flip b))(i (iden (length (car b))))(j 0))
     (if (and (null? k) (null? i)) '()
        (if (and (not (null? a)) (= j (car a)))
           (cons (car k) (puff (cdr a) (cdr k)      i  (+ j 1)))
           (cons (car i) (puff      a       k  (cdr i) (+ j 1)))
        ))))))
we can do inter thus:
(define (inter a b) (oss (rref (append (oss (rref a)) (oss (rref b))))))
I have copied the necessary Scheme fragments here to perform these calculations. Some code to test intersection code. See also the version in the repository. There is a demo at the end.

Unruly Hair

While (S⊥)⊥ = S indeed, there are subtleties. When the full vector space X and its dual X* are firmly in view and S is a subspace of X, then (S⊥)⊥ = S but what is the routine oss to do when it receives the 0 by n matrix for which to find the orthogonal space. A 0 by 3 matrix is identical to a 0 by 4 matrix. It the argument were known to be the 0 by 3 matrix then the correct response would be the 3 by 3 identity, which is not the same as the 4 by 4 identity. We could set the number of dimensions of X by parameterization but that is awkward and the value is seldom unavailable in the input values. Computing (X⊥)⊥ is thus a problem for (X⊥) is {0}. It is like expecting 1./(1./x) to yield the same as x even when x is zero.