QR decomposition

QR.tq

let qr mA =
  let m, n = Matrix.dimensions mA in
  let mI = Matrix.identity m in
  for 0 1 (n-1) (mI, mI, mA) [(mQ, mR, mQA), k ->
    let x = Matrix.column 0 mQA in
    let xb = Vector.norm x in
    let u = Array.mapi [i, xi -> if i=0 then xi - xb else xi] x in
    let v = Vector.scale (1 / Vector.norm u) u in
    let q k (i, j) =
      let x = if i=j then 1 else 0 in
      if i<k || j<k then x else
        x - 2 * Vector.get v (i-k) * Vector.get v (j-k) in
    let mQA = Matrix.mul (Matrix.init (m-k) (m-k) (q 0)) mQA in
    let mQA = Matrix.init (m-k-1) (n-k-1) [i, j ->
      Matrix.get mQA (i+1, j+1)] in
    let mQk = Matrix.init m m (q k) in
    Matrix.mul mQ (Matrix.transpose mQk), Matrix.mul mQk mR, mQA]
  @ [mQ, mR, _ -> mQ, Matrix.mul mR mA]