Intersection of Subspaces

We compute the intersection of two subspaces when subspaces are expressed by rows of spanning vectors. The dual space idea allows a concept of orthogonal space without an inner product.

We use Dirac’s notation <p|q> where q belongs to X and p belongs to X*, the dual space of X; this should keep the essential symmetry in view. <p|q> = pq. Recall that p is an 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 → <q|x> = 0} is a linear subspace of X*. (X^⊥)^⊥ = X and for subspaces U and V of X we have (U∩V) =(U^⊥ ∪ V^⊥)^⊥.

If {b_i} is a basis for X then there is a unique basis {c_i} for X* such that <c_j|b_i> = δ_ij.

If we have a subspace U of X spanned by a set u_i 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) =(U^⊥ ∪ V^⊥)^⊥ to compute the intersection of subspaces from their spanning sets.

Computing a spanning set for U^⊥

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 uniquely determines R. The row vectors of R also span U. If our original vectors were independent there will be no zero rows in R. Some of the columns of R will be all 0 except for one 1 and that 1 will be the first 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* is also in reduced row echelon form except the rows are reversed from top to bottom and the columns are reversed from right to left.
Column j of R is a sentinel if and only if column j of R* is not.
Extract the numbers from non-sentinel columns of R into a matrix, transpose that matrix, change the signs of the numbers; the columns of the new matrix are the non-sentinel columns of R*, in order.

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 <q|p> = 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; (dssp 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 (dssp 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) =(U^⊥ ∪ V^⊥)^⊥ we try the following:

(define (inter a b) (let* ((x (dssp (rref (minspan a))))(y (dssp (rref (minspan b))))
    (z (minspan (append x y)))) (dssp (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 dssp to ignore trailing zero vectors in its argument. With this version dss of dssp:

(define (dss 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) (dss (rref (append (dss (rref a)) (dss (rref b))))))

I have copied the necessary Scheme fragments here to perform these calculations. Some code to test intersection code. There is a demo at the end.

Computing a spanning set for U⊥

Computing a spanning set for U^⊥