QR decomposition

QR.clj

(require '[clojure.core.matrix :as m])

(defn mup [mat rs cs f]
  (m/set-selection mat rs cs
    (f (m/select mat rs cs))))

(defn qr
  "
    QR decomposition by Householder reflection
  
    adapted from Matlab implementation by @tobydriscoll

    returns Q the eigenvectors and R
    the diagonal of which is the eigenvalues

    Qd = I 0
         0 F

    F begins at column d, row d

    F is a n-d dimensional vector space
    in which a hyperplane H reflects z
    to the vector |z|e1

    v = |z|e1 - z

    Fy = (I - 2(vv'/v'v))y

    which is an orthonormal projector
  "
  ([a]
   (qr (m/identity-matrix (first (m/shape a))) (m/shape a) a))
  ([I [m n] A]
    (loop [d 0 R A Q I]
      (if (< d n)
        (let [
              [z1 :as z] (m/select R (range d m) d)
              v (m/matrix
                  (cons
                    (- (* -1.0 (Math/signum (double z1)) (m/magnitude z)) z1)
                    (m/mul -1.0 (m/select z :rest))))
              Qd (mup I (range d m) (range d n)
                    (fn Fy [i] (m/sub i (m/mul 2.0 (m/div (m/outer-product v v)
                                                          (m/inner-product v v))))))
              ]
          (recur (inc d) (m/mmul Qd R) (m/mmul Q Qd)))
        {:A A
         :Q Q
         :R R
         :A=QR (m/mmul Q R)}))))