Basis for the intersection of vector spaces

Question

The MuPAD Notebook Interface provides the linalg::intBasis function: http://www.mathworks.com/help/symbolic/mupad_ref/linalg-intbasis.html

How can I get the same functionality in Mathematica?

There exists an old thread called "Intersection of two vector spaces": Intersection of two vector spaces However, I am searching for a general solution which works with any number of vector spaces (like the linalg::intBasis does), not just two.

Moreover, the solution given in this old thread misbehaves in many cases (that I tried).

For three vector spaces I tried getIntersection[V1, getIntersection[V2, V3]] but it returned an incorrect result (where V1, V2 and V3 come from the "Example 1" from the linalg::intBasis web link given above).

Also, if I try getIntersection[V1, V1] it happily returns the {0, 0} vector among the basis vectors.

Additionally, the first comment (by Daniel Lichtblau) in the first "Answer" in this old thread says that one can also use coeffs=ker[[All,1;;n]]; ker.coeffs but the ker.coeffs part often generates an error: Dot::dotsh: Tensors {...} and {...} have incompatible shapes.

(Note: I'm a newcomer here and apparently I have no way to post "comments" in threads which do not belong to me as I get an error saying "You must have 50 reputation to comment" and my "answers" there get deleted. So, I am unable to report problems which I find in them.)

Update (2019.08.07): I have been using the code provided here for several years now. I have found that in some cases, especially when dealing with inexact numerical matrices, it could misbehave. I finally decided to do something about it and, with the help provided in this thread, I came out with this improved function, which may be useful for a casual trespasser who reads this. Note that the Chop function uses a default tolerance of 10.^(-10) but, you can easily change it in the getIntersectionBasis below, e.g.: use Chop[FullSimplify[#], 1.*^-14] (for a tolerance of 10.^(-14)).

getIntersectionBasis[] := {}
getIntersectionBasis[{}] := {}
getIntersectionBasis[{}, __] := {}
getIntersectionBasis[__, {}] := {}
getIntersectionBasis[l1_] := getIntersectionBasis[l1, l1]
getIntersectionBasis[l1_, l2_, l3__] := 
  getIntersectionBasis[l1, getIntersectionBasis[l2, l3]]
getIntersectionBasis[l1_, l2_] := 
  Catch[With[{ker = FullSimplify[NullSpace[FullSimplify[Transpose[Join[l1, l2]]], ZeroTest->(PossibleZeroQ[Chop[FullSimplify[#]]]&)]]}, 
    If[ker === {}, Throw[{}], 
      DeleteCases[FullSimplify[RowReduce[FullSimplify[ker[[All, 1 ;; Length[l1]]].l1], ZeroTest->(PossibleZeroQ[Chop[FullSimplify[#]]]&)]], {__?(PossibleZeroQ[Chop[FullSimplify[#]]]&)}]]] 
  ]

Answer 1

You'd need to row reduce the intersection set to remove linear dependencies it might have. The code below should handle this.

getIntersectionBasis[] := {}
getIntersectionBasis[{}] := {}
getIntersectionBasis[{}, __] := {}
getIntersectionBasis[__, {}] := {}
getIntersectionBasis[l1_] := getIntersectionBasis[l1, l1]
getIntersectionBasis[l1_, l2_, l3__] := 
  getIntersectionBasis[getIntersectionBasis[l1, l2], l3]
getIntersectionBasis[l1_, l2_] := 
  Catch[With[{ker = NullSpace[Transpose[Join[l1, l2]]]}, 
    If[ker === {}, Throw[{}], 
      DeleteCases[RowReduce[ker[[All, 1 ;; Length[l1]]].l1], {0 ..}]]]]

Now fill in your favorite example (which actually belongs in the question).