The Kitchen Sink and Other Oddities

Atabey Kaygun

Fast Null-Space Calculation via LU-Decomposition

Description of the Problem

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.

Some Theory

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.

An Implementation in Clojure

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/where

Let 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-space

Let 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]