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.
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:
( 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.