Today, I am going to work on a linear algebra problem. Let us assume we have a matrix \(A\), and I would like to calculate the null space of \(A\) (kernel of \(A\) for those mathematicians in the audience). I found this neat trick from MathOverflow with the MATLAB code that goes with it.
The LU decomposition splits a matrix as \[ A = PLU \] where \(L\) is an invertible lower triangular square matrix, \(U\) is an upper triangular matrix and \(P\) is a permutation matrix. Thus the rank of \(U\) determines the dimension of the null space, and the image \(L^{-1}P^{-1}\) determines the null space.
I am going to use JBlas to implement all the linear algebraic operations I am going to use today. First, let us load it up.
(defn where [pred coll]
(->> (map vector coll (range))
(filter #(pred (first %)))
(map second)))#'user/whereLet us test
(->> [4 4 0 0 3 3 1 1]
(where #(> % 2))
(into []))[0 1 4 5]Next, my implementation of the kernel of a matrix.
(defn null-space [mat err]
(let [LU (Decompose/lu mat)
JJ (->> LU .u MatrixFunctions/abs .rowMaxs .toArray seq
(where #(< % err)))
LL (.mmul (->> LU .l Solve/pinv)
(->> LU .p Solve/pinv))]
(->> JJ (map #(.getRow LL %)))))#'user/null-spaceLet me test:
(let [D (FloatMatrix/ones 4 4)]
(into [] (map #(.mmul % D) (null-space D 1e-4))))[#object[org.jblas.FloatMatrix 0x4c27d39d "[-0.000000, -0.000000, -0.000000, -0.000000]"] #object[org.jblas.FloatMatrix 0x40ee0a22 "[-0.000000, -0.000000, -0.000000, -0.000000]"] #object[org.jblas.FloatMatrix 0x7bde1f3a "[-0.000000, -0.000000, -0.000000, -0.000000]"]]Let me test on a larger example
(let [N (* (+ (rand-int 20) 1) 3)
D (take (/ (* N N) 3) (repeatedly rand))
X (as-> (concat D D D) $
(into-array Double/TYPE $)
(DoubleMatrix. $)
(.reshape $ N N))]
(into [] (map #(-> % (.mmul X) .norm2) (null-space X 1e-4))))[5.432236585530811E-15 5.979408099313226E-15 7.94990097832517E-15 5.952266062505376E-15 1.1430339784995241E-14 4.6055623193966644E-15 7.299039800550875E-15 8.471053469018959E-15 1.940347925441233E-14 3.473239481119178E-15 4.643016737199561E-15 4.1245726924716544E-15 9.721153472122608E-15 6.569709199941931E-15 2.183287502309909E-14 7.334201458919592E-15 2.6252268381282792E-14 1.541078123732933E-14 5.9509984362189764E-15 1.9378958691590616E-14 6.723954251254055E-15 8.887717976681934E-15 1.25297637180011E-14 6.841227998991561E-15 6.458786624653354E-15 7.171940574863952E-15 6.821762710232571E-15 7.565785875886503E-15]