class Matrix
Attributes
Public Class Methods
Creates a matrix n
x m
If you provide a block, it will be used to set the values. If not, val
will be used
# File lib/extendmatrix.rb, line 320 def build(n,m,val=0) f = (block_given?)? lambda {|i,j| yield(i, j)} : lambda {|i,j| val} Matrix.rows((0...n).collect {|i| (0...m).collect {|j| f.call(i,j)}}) end
Creates a matrix with the given matrices as diagonal blocks
# File lib/extendmatrix.rb, line 328 def diag(*args) dsize = 0 sizes = args.collect{|e| x = (e.is_a?(Matrix)) ? e.row_size : 1; dsize += x; x} m = Matrix.zero(dsize) count = 0 sizes.size.times{|i| range = count..(count+sizes[i]-1) m[range, range] = args[i] count += sizes[i] } m end
Tests if all the diagonal elements of two matrix are equal in delta
# File lib/extendmatrix.rb, line 359 def diag_in_delta?(m1, m0, eps = 1.0e-10) n = m1.row_size return false if n != m0.row_size or m1.column_size != m0.column_size n.times{|i| return false if (m1[i,i]-m0[i,i]).abs > eps } true end
Tests if all the elements of two matrix are equal in delta
# File lib/extendmatrix.rb, line 345 def equal_in_delta?(m0, m1, delta = 1.0e-10) delta = delta.abs m0.row_size.times {|i| m0.column_size.times {|j| x=m0[i,j]; y=m1[i,j] return false if (x < y - delta or x > y + delta) } } true end
Return an empty matrix on n=0 else, returns I(n)
# File lib/extendmatrix.rb, line 226 def self.robust_I(n) n>0 ? self.I(n) : self.[] end
Public Instance Methods
Return a value or a vector/matrix of values depending if the indexes are ranges or not m = Matrix.new_with_value(4, 3){|i, j| i * 3 + j} m: 0 1 2
3 4 5 6 7 8 9 10 11
m[1, 2] => 5 m => Vector[10, 11] m[0..1, 0..2] => Matrix[[0, 1, 2], [3, 4, 5]]
# File lib/extendmatrix.rb, line 241 def [](i, j) case i when Range case j when Range Matrix[*i.collect{|l| self.row(l)[j].to_a}] else column(j)[i] end else case j when Range row(i)[j] else ids(i, j) end end end
Set the values of a matrix m = Matrix.build
(3, 3){|i, j| i * 3 + j} m: 0 1 2
3 4 5 6 7 8
m[1, 2] = 9 => Matrix[[0, 1, 2], [3, 4, 9], [6, 7, 8]] m = Vector[8, 8] => Matrix[[0, 1, 2], [3, 8, 8], [6, 7, 8]] m[0..1, 0..1] = Matrix[[0, 0, 0],[0, 0, 0]]
=> Matrix[[0, 0, 2], [0, 0, 8], [6, 7, 8]]
# File lib/extendmatrix.rb, line 270 def []=(i, j, v) case i when Range if i.entries.size == 1 self[i.begin, j] = (v.is_a?(Matrix) ? v.row(0) : v) else case j when Range if j.entries.size == 1 self[i, j.begin] = (v.is_a?(Matrix) ? v.column(0) : v) else i.each{|l| self.row= l, v.row(l - i.begin), j} end else self.column= j, v, i end end else case j when Range if j.entries.size == 1 self[i, j.begin] = (v.is_a?(Vector) ? v[0] : v) else self.row= i, v, j end else @rows[i][j] = (v.is_a?(Vector) ? v[0] : v) end end end
Returns a list with the maximum lengths
# File lib/extendmatrix.rb, line 429 def cols_len (0...column_size).collect {|j| max_len_column(j)} end
Return a duplicate matrix, with all elements copied
# File lib/extendmatrix.rb, line 304 def dup (self.class).rows(self.rows.dup) end
Iterate the elements of a matrix
# File lib/extendmatrix.rb, line 463 def each @rows.each {|x| x.each {|e| yield(e)}} nil end
# File lib/extendmatrix.rb, line 308 def initialize_copy(orig) init_rows(orig.rows, true) self.wrap=(orig.wrap) end
Returns the maximum length of column elements
# File lib/extendmatrix.rb, line 422 def max_len_column(j) column_collect(j) {|x| x.to_s.length}.max end
Returns the matrix divided by a scalar
# File lib/extendmatrix.rb, line 372 def quo(v) map {|e| e.quo(v)} end
Set de values of a matrix and the value of wrap property
# File lib/extendmatrix.rb, line 385 def set(m) 0.upto(m.row_size - 1) do |i| 0.upto(m.column_size - 1) do |j| self[i, j] = m[i, j] end end self.wrap = m.wrap end
Returns a string for nice printing matrix
# File lib/extendmatrix.rb, line 437 def to_s(mode = :pretty, len_col = 3) return "Matrix[]" if empty? if mode == :pretty clen = cols_len to_a.collect {|r| i=0 r.map {|x| l=clen[i] i+=1 format("%#{l}s ", x.to_s) } << "\n" }.join("") else i = 0; s = ""; cs = column_size each do |e| i = (i + 1) % cs s += format("%#{len_col}s ", e.to_s) s += "\n" if i == 0 end s end end
Set wrap feature of a matrix
# File lib/extendmatrix.rb, line 409 def wrap=(mode = :torus) case mode when :torus then eval(wraplate("i %= row_size; j %= column_size")) when :h_cylinder then eval(wraplate("i %= row_size")) when :v_cylinder then eval(wraplate("j %= column_size")) when :nil then eval(wraplate) end @wrap = mode end
# File lib/extendmatrix.rb, line 394 def wraplate(ijwrap = "") "class << self def [](i, j) #{ijwrap}; @rows[i][j] end def []=(i, j, v) #{ijwrap}; @rows[i][j] = v end end" end
Advanced methods
↑ topConstants
- EXTENSION_VERSION
Public Instance Methods
Returns the upper bidiagonal matrix obtained with Householder
Bidiagonalization algorithm
# File lib/extendmatrix.rb, line 825 def bidiagonal ub, vb = Householder.bidiag(self) ub.t * self * vb end
Classical Jacobi
8.4.3 Golub & van Loan
# File lib/extendmatrix.rb, line 1074 def cJacobi(tol = 1.0e-10) a = self.clone n = row_size v = Matrix.I(n) eps = tol * a.normF while Jacobi.off(a) > eps p, q = Jacobi.max(a) c, s = Jacobi.sym_schur2(a, p, q) #print "\np:#{p} q:#{q} c:#{c} s:#{s}\n" j = Jacobi.J(p, q, c, s, n) a = j.t * a * j v = v * j end return a, v end
Returns the aproximation matrix computed with Classical Jacobi
algorithm. The aproximate eigenvalues values are in the diagonal of the matrix A.
# File lib/extendmatrix.rb, line 1094 def cJacobiA(tol = 1.0e-10) cJacobi(tol)[0] end
Returns the orthogonal matrix obtained with the Jacobi
eigenvalue algorithm. The columns of V are the eigenvector.
# File lib/extendmatrix.rb, line 1111 def cJacobiV(tol = 1.0e-10) cJacobi(tol)[1] end
Returns column vector number “j” as a Vector
. When the block is given, the elements of column “j” are mmodified
# File lib/extendmatrix.rb, line 591 def column!(j) if block_given? (0...row_size).collect { |i| @rows[i][j] = yield @rows[i][j] } else column(j) end end
Returns the column/s of matrix as a Matrix
# File lib/extendmatrix.rb, line 665 def column2matrix(c) if c.is_a?(Range) a=c.map {|v| column(v).to_a}.find_all {|v| v.size>0} else a = column(c).to_a end if c.is_a?(Range) return Matrix.columns(a) else return Matrix[*a.collect{|x| [x]}] end end
Set a certain column with the values of a Vector
m = Matrix.build
(3, 3){|i, j| i * 3 + j + 1} m.column= 1, Vector[1, 1, 1], 1..2 m => 1 2 3
4 1 6 7 1 9
# File lib/extendmatrix.rb, line 608 def column=(args) m = row_size c, v, r = MMatrix.id_vect_range(args, m) (m..r.begin - 1).each{|i| self[i, c] = 0} [v.size, r.entries.size].min.times{|i| self[i + r.begin, c] = v[i]} ((v.size + r.begin)..r.entries.last).each {|i| self[i, c] = 0} end
Returns an array with the elements collected from the column “j”. When a block is given, the elements of that vector are iterated.
# File lib/extendmatrix.rb, line 582 def column_collect(j, &block) f = MMatrix.default_block(block) (0...row_size).collect {|r| f.call(self[r, j])} end
Returns the marginal of columns
# File lib/extendmatrix.rb, line 685 def column_sum (0...column_size).collect {|i| column(i).sum } end
Tests if the matrix is empty or not
# File lib/extendmatrix.rb, line 641 def empty? @rows.empty? if @rows end
Returns the upper triunghiular matrix R of a Givens
QR factorization
# File lib/extendmatrix.rb, line 932 def givensR Givens.QR(self)[0] end
Modified Gram Schmidt QR factorization (MC, Golub, p. 232) A = Q_1 * R_1
# File lib/extendmatrix.rb, line 856 def gram_schmidt a = clone n = column_size m = row_size q = Matrix.build(m, n){0} r = Matrix.zero(n) for k in 0...n r[k,k] = a[0...m, k].norm q[0...m, k] = a[0...m, k] / r[k, k] for j in (k+1)...n r[k, j] = q[0...m, k].t * a[0...m, j] a[0...m, j] -= q[0...m, k] * r[k, j] end end return q, r end
Returns the Q_1 matrix of Modified Gram Schmidt algorithm Q_1 has orthonormal columns
# File lib/extendmatrix.rb, line 877 def gram_schmidtQ gram_schmidt[0] end
Returns the R_1 upper triangular matrix of Modified Gram Schmidt algorithm
# File lib/extendmatrix.rb, line 884 def gram_schmidtR gram_schmidt[1] end
Returns the orthogonal matrix Q of Hessenberg
QR factorization Q = G_1 … G_(n-1) where G_j is the Givens
rotation G_j = G(j, j+1, omega_j)
# File lib/extendmatrix.rb, line 967 def hessenbergQ Hessenberg.QR(self)[0] end
Returns the upper triunghiular matrix R of a Hessenberg
QR factorization
# File lib/extendmatrix.rb, line 974 def hessenbergR Hessenberg.QR(self)[1] end
Return an upper Hessenberg
matrix obtained with Householder
reduction to Hessenberg
Form algorithm
# File lib/extendmatrix.rb, line 981 def hessenberg_form_H Householder.toHessenberg(self)[0] end
Returns the orthogonal matrix Q of Householder
QR factorization where Q = H_1 * H_2 * H_3 * … * H_n,
# File lib/extendmatrix.rb, line 834 def houseQ h = Householder.QR(self) q = h[0] (1...h.size).each{|i| q *= h[i]} q end
Returns the matrix R of Householder
QR factorization R = H_n * H_n-1 * … * H_1 * A is an upper triangular matrix
# File lib/extendmatrix.rb, line 845 def houseR h = Householder.QR(self) r = self.clone h.size.times{|i| r = h[i] * r} r end
# File lib/extendmatrix.rb, line 629 def norm(p = 2) Vector::Norm.sqnorm(self, p) ** (1.quo(p)) end
# File lib/extendmatrix.rb, line 633 def norm_frobenius norm end
The real Schur decomposition. The eigenvalues are aproximated in diagonal elements of the real Schur decomposition matrix
# File lib/extendmatrix.rb, line 989 def realSchur(eps = 1.0e-10, steps = 100) h = self.hessenberg_form_H h1 = Matrix[] i = 0 loop do h1 = h.hessenbergR * h.hessenbergQ break if Matrix.diag_in_delta?(h1, h, eps) or steps <= 0 h = h1.clone steps -= 1 i += 1 end h1 end
Returns row vector number “i” like Matrix.row as a Vector
. When the block is given, the elements of row “i” are modified
# File lib/extendmatrix.rb, line 570 def row!(i) if block_given? @rows[i].collect! {|e| yield e } else Vector.elements(@rows[i], false) end end
Returns row(s) of matrix as a Matrix
# File lib/extendmatrix.rb, line 649 def row2matrix(r) if r.is_a? Range a=r.map {|v| v<row_size ? row(v).to_a : nil }.find_all {|v| !v.nil?} else a = row(r).to_a end if r.is_a?(Range) and r.entries.size > 1 return Matrix[*a] else return Matrix[a] end end
Set a certain row with the values of a Vector
m = Matrix.new(3, 3){|i, j| i * 3 + j + 1} m.row= 0, Vector[0, 0], 1..2 m => 1 0 0
4 5 6 7 8 9
# File lib/extendmatrix.rb, line 624 def row=(args) i, val, range = MMatrix.id_vect_range(args, column_size) row!(i)[range] = val end
Returns an array with the elements collected from the row “i”. When a block is given, the elements of that vector are iterated.
# File lib/extendmatrix.rb, line 562 def row_collect(i, &block) f = MMatrix.default_block(block) @rows[i].collect {|e| f.call(e)} end
Returns the marginal of rows
# File lib/extendmatrix.rb, line 679 def row_sum (0...row_size).collect {|i| row(i).sum } end
Calculate sum of cells
# File lib/extendmatrix.rb, line 691 def total_sum row_sum.inject(&:+) end
SPSS like methods
↑ topPublic Instance Methods
Diagonal of matrix. Returns an array with as many elements as the minimum of the number of rows and the number of columns in the argument
# File lib/extendmatrix.rb, line 526 def diagonal min = (row_size<column_size) ? row_size : column_size min.times.map {|i| self[i,i]} end
Element wise multiplication
# File lib/extendmatrix.rb, line 479 def e_mult(other) elementwise_operation(:*,other) end
Element wise multiplication
# File lib/extendmatrix.rb, line 484 def e_quo(other) elementwise_operation(:quo,other) end
Returns eigenvalues and eigenvectors of a matrix on a Hash like CALL EIGEN on SPSS.
-
:eigenvectors: contains the eigenvectors as columns of a new
Matrix
, ordered in descendent order -
:eigenvalues: contains an array with the eigenvalues, ordered in descendent order
# File lib/extendmatrix.rb, line 505 def eigen ep=eigenpairs { :eigenvalues=>ep.map {|a| a[0]} , :eigenvectors=>Matrix.columns(ep.map{|a| a[1]}) } end
# File lib/extendmatrix.rb, line 493 def eigenpairs eigval, eigvec= eigenvaluesJacobi, cJacobiV eigenpairs=eigval.size.times.map {|i| [eigval[i], eigvec.column(i)] } eigenpairs=eigenpairs.sort{|a,b| a[0]<=>b[0]}.reverse end
Element wise operation
# File lib/extendmatrix.rb, line 472 def elementwise_operation(op,other) Matrix.build(row_size,column_size) do |row, column| self[row,column].send(op,other[row,column]) end end
Returns a new matrix with the natural logarithm of initial values
# File lib/extendmatrix.rb, line 512 def ln map {|v| Math::log(v)} end
Matrix
sum of squares
# File lib/extendmatrix.rb, line 489 def mssq @rows.inject(0){|ac,row| ac+(row.inject(0) {|acr,i| acr+(i**2)})} end
Returns a new matrix, with the square root of initial values
# File lib/extendmatrix.rb, line 516 def sqrt map {|v| Math::sqrt(v)} end
Sum of squares and cross-products. Equal to m.t * m
# File lib/extendmatrix.rb, line 520 def sscp transpose*self end