Class: Matrix

Inherits:
Object
  • Object
show all
Includes:
Enumerable
Defined in:
lib/version.rb,
lib/extendmatrix.rb

Defined Under Namespace

Modules: Givens, Hessenberg, Householder, Jacobi, LU, MMatrix

Constant Summary collapse

EXTENSION_VERSION = 
"0.4"

Instance Attribute Summary collapse

Class Method Summary collapse

Instance Method Summary collapse

Instance Attribute Details

#rows ⇒ Object (readonly)

Returns the value of attribute rows.



216
217
218
 
# File 'lib/extendmatrix.rb', line 216

def rows
  @rows
end

#wrap ⇒ Object

Returns the value of attribute wrap.



216
217
218
 
# File 'lib/extendmatrix.rb', line 216

def wrap
  @wrap
end

Class Method Details

.diag(*args) ⇒ Object

Creates a matrix with the given matrices as diagonal blocks



328
329
330
331
332
333
334
335
336
337
338
339
340
 
# 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

.diag_in_delta?(m1, m0, eps = 1.0e-10) ⇒ Boolean

Tests if all the diagonal elements of two matrix are equal in delta

Returns:

  • (Boolean) 


359
360
361
362
363
364
365
366
 
# 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

.equal_in_delta?(m0, m1, delta = 1.0e-10) ⇒ Boolean

Tests if all the elements of two matrix are equal in delta

Returns:

  • (Boolean) 


345
346
347
348
349
350
351
352
353
354
 
# 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

.robust_I(n) ⇒ Object

Return an empty matrix on n=0 else, returns I(n)



226
227
228
 
# File 'lib/extendmatrix.rb', line 226

def self.robust_I(n)
  n>0 ? self.I(n) : self.[]
end

Instance Method Details

#[](i, j) ⇒ Object

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]]



241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
 
# 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

#[]=(i, j, v) ⇒ Object

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]]



270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
 
# 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

#bidiagonal ⇒ Object

Returns the upper bidiagonal matrix obtained with Householder Bidiagonalization algorithm



825
826
827
828
 
# File 'lib/extendmatrix.rb', line 825

def bidiagonal
  ub, vb = Householder.bidiag(self)
  ub.t * self * vb
end

#build(n, m, val = 0) ⇒ Object

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



320
321
322
323
 
# 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

#cJacobi(tol = 1.0e-10) ⇒ Object

Classical Jacobi 8.4.3 Golub & van Loan



1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
 
# 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

#cJacobiA(tol = 1.0e-10) ⇒ Object

Returns the aproximation matrix computed with Classical Jacobi algorithm. The aproximate eigenvalues values are in the diagonal of the matrix A.



1094
1095
1096
 
# File 'lib/extendmatrix.rb', line 1094

def cJacobiA(tol = 1.0e-10)
  cJacobi(tol)[0]
end

#cJacobiV(tol = 1.0e-10) ⇒ Object

Returns the orthogonal matrix obtained with the Jacobi eigenvalue algorithm. The columns of V are the eigenvector.



1111
1112
1113
 
# File 'lib/extendmatrix.rb', line 1111

def cJacobiV(tol = 1.0e-10)
  cJacobi(tol)[1]
end

#cols_len ⇒ Object

Returns a list with the maximum lengths



429
430
431
 
# File 'lib/extendmatrix.rb', line 429

def cols_len
  (0...column_size).collect {|j| max_len_column(j)}
end

#column!(j) ⇒ Object Also known as: column_collect!

Returns column vector number “j” as a Vector. When the block is given, the elements of column “j” are mmodified



591
592
593
594
595
596
597
 
# 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

#column2matrix(c) ⇒ Object

Returns the column/s of matrix as a Matrix



665
666
667
668
669
670
671
672
673
674
675
676
 
# 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

#column=(args) ⇒ Object

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



608
609
610
611
612
613
614
 
# 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

#column_collect(j, &block) ⇒ Object

Returns an array with the elements collected from the column “j”. When a block is given, the elements of that vector are iterated.



582
583
584
585
 
# 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

#column_sum ⇒ Object

Returns the marginal of columns



685
686
687
688
689
 
# File 'lib/extendmatrix.rb', line 685

def column_sum
  (0...column_size).collect {|i|
    column(i).sum
  }
end

#diagonal ⇒ Object

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



526
527
528
529
 
# 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

#dup ⇒ Object

Return a duplicate matrix, with all elements copied



304
305
306
 
# File 'lib/extendmatrix.rb', line 304

def dup
  (self.class).rows(self.rows.dup)
end

#e_mult(other) ⇒ Object

Element wise multiplication



479
480
481
 
# File 'lib/extendmatrix.rb', line 479

def e_mult(other)
  elementwise_operation(:*,other)
end

#e_quo(other) ⇒ Object

Element wise multiplication



484
485
486
 
# File 'lib/extendmatrix.rb', line 484

def e_quo(other)
  elementwise_operation(:quo,other)
end

#each ⇒ Object

Iterate the elements of a matrix



463
464
465
466
 
# File 'lib/extendmatrix.rb', line 463

def each
  @rows.each {|x| x.each {|e| yield(e)}}
  nil
end

#eigen ⇒ Object

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



505
506
507
508
509
510
 
# 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

#eigenpairs ⇒ Object 



493
494
495
496
497
498
499
 
# 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

#eigenvaluesJacobi ⇒ Object

Returns a Vector with the eigenvalues aproximated values. The eigenvalues are computed with the Classic Jacobi Algorithm.



1102
1103
1104
1105
 
# File 'lib/extendmatrix.rb', line 1102

def eigenvaluesJacobi
  a = cJacobiA
  Vector[*(0...row_size).collect{|i| a[i, i]}]
end

#elementwise_operation(op, other) ⇒ Object

Element wise operation



472
473
474
475
476
 
# 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

#empty? ⇒ Boolean

Tests if the matrix is empty or not

Returns:

  • (Boolean) 


641
642
643
 
# File 'lib/extendmatrix.rb', line 641

def empty?
  @rows.empty? if @rows
end

#givensQ ⇒ Object

Returns the orthogonal matrix Q of Givens QR factorization. Q = G_1 * … * G_t where G_j is the j’th Givens rotation and ‘t’ is the total number of rotations



941
942
943
 
# File 'lib/extendmatrix.rb', line 941

def givensQ
  Givens.QR(self)[1]
end

#givensR ⇒ Object

Returns the upper triunghiular matrix R of a Givens QR factorization



932
933
934
 
# File 'lib/extendmatrix.rb', line 932

def givensR
  Givens.QR(self)[0]
end

#gram_schmidt ⇒ Object

Modified Gram Schmidt QR factorization (MC, Golub, p. 232) A = Q_1 * R_1



856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
 
# 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

#gram_schmidtQ ⇒ Object

Returns the Q_1 matrix of Modified Gram Schmidt algorithm Q_1 has orthonormal columns



877
878
879
 
# File 'lib/extendmatrix.rb', line 877

def gram_schmidtQ
  gram_schmidt[0]
end

#gram_schmidtR ⇒ Object

Returns the R_1 upper triangular matrix of Modified Gram Schmidt algorithm



884
885
886
 
# File 'lib/extendmatrix.rb', line 884

def gram_schmidtR
  gram_schmidt[1]
end

#hessenberg_form_H ⇒ Object

Return an upper Hessenberg matrix obtained with Householder reduction to Hessenberg Form algorithm



981
982
983
 
# File 'lib/extendmatrix.rb', line 981

def hessenberg_form_H
  Householder.toHessenberg(self)[0]
end

#hessenbergQ ⇒ Object

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)



967
968
969
 
# File 'lib/extendmatrix.rb', line 967

def hessenbergQ
  Hessenberg.QR(self)[0]
end

#hessenbergR ⇒ Object

Returns the upper triunghiular matrix R of a Hessenberg QR factorization



974
975
976
 
# File 'lib/extendmatrix.rb', line 974

def hessenbergR
  Hessenberg.QR(self)[1]
end

#houseQ ⇒ Object

Returns the orthogonal matrix Q of Householder QR factorization where Q = H_1 * H_2 * H_3 * … * H_n,



834
835
836
837
838
839
 
# 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

#houseR ⇒ Object

Returns the matrix R of Householder QR factorization R = H_n * H_n-1 * … * H_1 * A is an upper triangular matrix



845
846
847
848
849
850
 
# 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

#ids ⇒ Object 



229
 
# File 'lib/extendmatrix.rb', line 229

alias :ids :[]

#initialize_copy(orig) ⇒ Object 



308
309
310
311
 
# File 'lib/extendmatrix.rb', line 308

def initialize_copy(orig)
  init_rows(orig.rows, true)
  self.wrap=(orig.wrap)
end

#initialize_old(init_method, *argv) ⇒ Object

For invoking a method



221
222
223
 
# File 'lib/extendmatrix.rb', line 221

def initialize_old(init_method, *argv) #:nodoc:
  self.send(init_method, *argv)
end

#L ⇒ Object

Return the lower triangular matrix of LU factorization L = M_1^-1 * … * M_n-1^-1 = I + sum_k=1^n-1 tau * e



744
745
746
 
# File 'lib/extendmatrix.rb', line 744

def L
  LU.factorization(self)[0]
end

#ln ⇒ Object

Returns a new matrix with the natural logarithm of initial values



512
513
514
 
# File 'lib/extendmatrix.rb', line 512

def ln
  map {|v| Math::log(v)}
end

#max_len_column(j) ⇒ Object

Returns the maximum length of column elements



422
423
424
 
# File 'lib/extendmatrix.rb', line 422

def max_len_column(j)
  column_collect(j) {|x| x.to_s.length}.max
end

#mssq ⇒ Object

Matrix sum of squares



489
490
491
 
# File 'lib/extendmatrix.rb', line 489

def mssq
  @rows.inject(0){|ac,row| ac+(row.inject(0) {|acr,i| acr+(i**2)})}
end

#norm(p = 2) ⇒ Object 



629
630
631
 
# File 'lib/extendmatrix.rb', line 629

def norm(p = 2)
  Vector::Norm.sqnorm(self, p) ** (1.quo(p))
end

#norm_frobenius ⇒ Object Also known as: normF 



633
634
635
 
# File 'lib/extendmatrix.rb', line 633

def norm_frobenius
  norm
end

#quo(v) ⇒ Object Also known as: /

Returns the matrix divided by a scalar



372
373
374
 
# File 'lib/extendmatrix.rb', line 372

def quo(v)
  map {|e| e.quo(v)}
end

#realSchur(eps = 1.0e-10, steps = 100) ⇒ Object

The real Schur decomposition. The eigenvalues are aproximated in diagonal elements of the real Schur decomposition matrix



989
990
991
992
993
994
995
996
997
998
999
1000
1001
 
# 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

#row!(i) ⇒ Object Also known as: row_collect!

Returns row vector number “i” like Matrix.row as a Vector. When the block is given, the elements of row “i” are modified



570
571
572
573
574
575
576
 
# 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

#row2matrix(r) ⇒ Object

Returns row(s) of matrix as a Matrix



649
650
651
652
653
654
655
656
657
658
659
660
 
# 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

#row=(args) ⇒ Object

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



624
625
626
627
 
# File 'lib/extendmatrix.rb', line 624

def row=(args)
  i, val, range = MMatrix.id_vect_range(args, column_size)
  row!(i)[range] = val
end

#row_collect(i, &block) ⇒ Object

Returns an array with the elements collected from the row “i”. When a block is given, the elements of that vector are iterated.



562
563
564
565
 
# File 'lib/extendmatrix.rb', line 562

def row_collect(i, &block)
  f = MMatrix.default_block(block)
  @rows[i].collect {|e| f.call(e)}
end

#row_sum ⇒ Object

Returns the marginal of rows



679
680
681
682
683
 
# File 'lib/extendmatrix.rb', line 679

def row_sum
  (0...row_size).collect {|i|
    row(i).sum
  }
end

#set(m) ⇒ Object

Set de values of a matrix and the value of wrap property



385
386
387
388
389
390
391
392
 
# 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

#sqrt ⇒ Object

Returns a new matrix, with the square root of initial values



516
517
518
 
# File 'lib/extendmatrix.rb', line 516

def sqrt
  map {|v| Math::sqrt(v)}
end

#sscp ⇒ Object

Sum of squares and cross-products. Equal to m.t * m



520
521
522
 
# File 'lib/extendmatrix.rb', line 520

def sscp
  transpose*self
end

#to_s(mode = :pretty, len_col = 3) ⇒ Object

Returns a string for nice printing matrix



437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
 
# 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

#to_s_old ⇒ Object 



432
 
# File 'lib/extendmatrix.rb', line 432

alias :to_s_old :to_s

#total_sum ⇒ Object

Calculate sum of cells



691
692
693
 
# File 'lib/extendmatrix.rb', line 691

def total_sum
  row_sum.inject(&:+)
end

#U ⇒ Object

Return the upper triangular matrix of LU factorization M_n-1 * … * M_1 * A = U



737
738
739
 
# File 'lib/extendmatrix.rb', line 737

def U
  LU.factorization(self)[1]
end

#wraplate(ijwrap = "") ⇒ Object 



394
395
396
397
398
399
400
401
402
403
404
 
# 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