Skip to content

Instantly share code, notes, and snippets.

@Hendekagon
Created September 12, 2021 09:09
Show Gist options
  • Save Hendekagon/08e1f1582df5a0423aefd70879d6ceff to your computer and use it in GitHub Desktop.
Save Hendekagon/08e1f1582df5a0423aefd70879d6ceff to your computer and use it in GitHub Desktop.
QR decomposition
(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)}))))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment