module Matrix::Householder
Public Class Methods
QR(mat)
click to toggle source
a QR
factorization that uses Householder
transformation Q^T * A = R MC, Golub & van Loan, pg 224, 5.2.1 Householder
QR
# File lib/extendmatrix.rb, line 754 def self.QR(mat) h = [] a = mat.clone m = a.row_size n = a.column_size n.times do |j| v, beta = a[j..m - 1, j].house h[j] = Matrix.diag(Matrix.robust_I(j), Matrix.I(m-j)- beta * (v * v.t)) a[j..m-1, j..n-1] = (Matrix.I(m-j) - beta * (v * v.t)) * a[j..m-1, j..n-1] a[(j+1)..m-1,j] = v[2..(m-j)] if j < m - 1 end h end
bidiag(mat)
click to toggle source
Householder
Bidiagonalization algorithm. MC, Golub, pg 252, Algorithm 5.4.2 Returns the matrices U_B and V_B such that: U_B^T * A * V_B = B, where B is upper bidiagonal.
# File lib/extendmatrix.rb, line 783 def self.bidiag(mat) a = mat.clone m = a.row_size n = a.column_size ub = Matrix.I(m) vb = Matrix.I(n) n.times do |j| v, beta = a[j..m-1,j].house a[j..m-1, j..n-1] = (Matrix.I(m-j) - beta * (v * v.t)) * a[j..m-1, j..n-1] a[j+1..m-1, j] = v[1..(m-j-1)] ub *= bidiagUV(a[j+1..m-1,j], m, beta) #Ub = U_1 * U_2 * ... * U_n if j < n - 2 v, beta = (a[j, j+1..n-1]).house a[j..m-1, j+1..n-1] = a[j..m-1, j+1..n-1] * (Matrix.I(n-j-1) - beta * (v * v.t)) a[j, j+2..n-1] = v[1..n-j-2] vb *= bidiagUV(a[j, j+2..n-1], n, beta) #Vb = V_1 * U_2 * ... * V_n-2 end end return ub, vb end
bidiagUV(essential, dim, beta)
click to toggle source
From the essential part of Householder
vector it returns the coresponding upper(U_j)/lower(V_j) matrix
# File lib/extendmatrix.rb, line 772 def self.bidiagUV(essential, dim, beta) v = Vector.concat(Vector[1], essential) dimv = v.size Matrix.diag(Matrix.robust_I(dim - dimv), Matrix.I(dimv) - beta * (v * v.t) ) end
toHessenberg(mat)
click to toggle source
Householder
Reduction to Hessenberg
Form
# File lib/extendmatrix.rb, line 807 def self.toHessenberg(mat) h = mat.clone n = h.row_size u0 = Matrix.I(n) for k in (0...n - 2) v, beta = h[k+1..n-1, k].house #the householder matrice part houseV = Matrix.I(n-k-1) - beta * (v * v.t) u0 *= Matrix.diag(Matrix.I(k+1), houseV) h[k+1..n-1, k..n-1] = houseV * h[k+1..n-1, k..n-1] h[0..n-1, k+1..n-1] = h[0..n-1, k+1..n-1] * houseV end return h, u0 end