Class: Matrix

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

Defined Under Namespace

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

Constant Summary collapse

EXTENSION_VERSION = 
"0.3.1"

Instance Attribute Summary collapse

Class Method Summary collapse

Instance Method Summary collapse

Instance Attribute Details

#rows ⇒ Object (readonly)

Returns the value of attribute rows.



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

def rows
  @rows
end

#wrap ⇒ Object

Returns the value of attribute wrap.



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

def wrap
  @wrap
end

Class Method Details

.diag(*args) ⇒ Object

Creates a matrix with the given matrices as diagonal blocks



348
349
350
351
352
353
354
355
356
357
358
359
360
 
# File 'lib/extendmatrix.rb', line 348

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) 


379
380
381
382
383
384
385
386
 
# File 'lib/extendmatrix.rb', line 379

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) 


365
366
367
368
369
370
371
372
373
374
 
# File 'lib/extendmatrix.rb', line 365

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)



238
239
240
 
# File 'lib/extendmatrix.rb', line 238

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



254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
 
# File 'lib/extendmatrix.rb', line 254

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



287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
 
# File 'lib/extendmatrix.rb', line 287

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



850
851
852
853
 
# File 'lib/extendmatrix.rb', line 850

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



340
341
342
343
 
# File 'lib/extendmatrix.rb', line 340

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



1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
 
# File 'lib/extendmatrix.rb', line 1095

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.



1115
1116
1117
 
# File 'lib/extendmatrix.rb', line 1115

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.



1132
1133
1134
 
# File 'lib/extendmatrix.rb', line 1132

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

#cols_len ⇒ Object

Returns a list with the maximum lengths



449
450
451
 
# File 'lib/extendmatrix.rb', line 449

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



612
613
614
615
616
617
618
 
# File 'lib/extendmatrix.rb', line 612

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



686
687
688
689
690
691
692
693
694
695
696
697
 
# File 'lib/extendmatrix.rb', line 686

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



629
630
631
632
633
634
635
 
# File 'lib/extendmatrix.rb', line 629

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.



603
604
605
606
 
# File 'lib/extendmatrix.rb', line 603

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



706
707
708
709
710
 
# File 'lib/extendmatrix.rb', line 706

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



542
543
544
545
 
# File 'lib/extendmatrix.rb', line 542

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



322
323
324
 
# File 'lib/extendmatrix.rb', line 322

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

#e_mult(other) ⇒ Object

Element wise multiplication



499
500
501
 
# File 'lib/extendmatrix.rb', line 499

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

#e_quo(other) ⇒ Object

Element wise multiplication



503
504
505
 
# File 'lib/extendmatrix.rb', line 503

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

#each ⇒ Object

Iterate the elements of a matrix



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

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



521
522
523
524
525
526
 
# File 'lib/extendmatrix.rb', line 521

def eigen
  ep=eigenpairs
  { :eigenvalues=>ep.map {|a| a[0]} ,
    :eigenvectors=>Matrix.columns(ep.map{|a| a[1]})
  }
end

#eigenpairs ⇒ Object 



510
511
512
513
514
515
516
 
# File 'lib/extendmatrix.rb', line 510

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.



1123
1124
1125
1126
 
# File 'lib/extendmatrix.rb', line 1123

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

#elementwise_operation(op, other) ⇒ Object

Element wise operation



493
494
495
496
497
 
# File 'lib/extendmatrix.rb', line 493

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) 


662
663
664
 
# File 'lib/extendmatrix.rb', line 662

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



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

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

#givensR ⇒ Object

Returns the upper triunghiular matrix R of a Givens QR factorization



957
958
959
 
# File 'lib/extendmatrix.rb', line 957

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

#gram_schmidt ⇒ Object

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



881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
 
# File 'lib/extendmatrix.rb', line 881

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



902
903
904
 
# File 'lib/extendmatrix.rb', line 902

def gram_schmidtQ
  gram_schmidt[0]
end

#gram_schmidtR ⇒ Object

Returns the R_1 upper triangular matrix of Modified Gram Schmidt algorithm



909
910
911
 
# File 'lib/extendmatrix.rb', line 909

def gram_schmidtR
  gram_schmidt[1]
end

#hessenberg_form_H ⇒ Object

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



1006
1007
1008
 
# File 'lib/extendmatrix.rb', line 1006

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)



992
993
994
 
# File 'lib/extendmatrix.rb', line 992

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

#hessenbergR ⇒ Object

Returns the upper triunghiular matrix R of a Hessenberg QR factorization



999
1000
1001
 
# File 'lib/extendmatrix.rb', line 999

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,



859
860
861
862
863
864
 
# File 'lib/extendmatrix.rb', line 859

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



870
871
872
873
874
875
 
# File 'lib/extendmatrix.rb', line 870

def houseR
  h = Householder.QR(self)
  r = self.clone
  h.size.times{|i| r = h[i] * r}
  r
end

#ids ⇒ Object 



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

alias :ids :[]

#initialize_copy(orig) ⇒ Object 



326
327
328
329
 
# File 'lib/extendmatrix.rb', line 326

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

#initialize_old(init_method, *argv) ⇒ Object

For invoking a method



233
234
235
 
# File 'lib/extendmatrix.rb', line 233

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



769
770
771
 
# File 'lib/extendmatrix.rb', line 769

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

#ln ⇒ Object

Returns a new matrix with the natural logarithm of initial values



528
529
530
 
# File 'lib/extendmatrix.rb', line 528

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

#max_len_column(j) ⇒ Object

Returns the maximum length of column elements



442
443
444
 
# File 'lib/extendmatrix.rb', line 442

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

#mssq ⇒ Object

Matrix sum of squares



507
508
509
 
# File 'lib/extendmatrix.rb', line 507

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

#norm(p = 2) ⇒ Object 



650
651
652
 
# File 'lib/extendmatrix.rb', line 650

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

#norm_frobenius ⇒ Object Also known as: normF 



654
655
656
 
# File 'lib/extendmatrix.rb', line 654

def norm_frobenius
  norm
end

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

Returns the matrix divided by a scalar



392
393
394
 
# File 'lib/extendmatrix.rb', line 392

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



1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
 
# File 'lib/extendmatrix.rb', line 1014

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



590
591
592
593
594
595
596
 
# File 'lib/extendmatrix.rb', line 590

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



670
671
672
673
674
675
676
677
678
679
680
681
 
# File 'lib/extendmatrix.rb', line 670

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



645
646
647
648
 
# File 'lib/extendmatrix.rb', line 645

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.



581
582
583
584
 
# File 'lib/extendmatrix.rb', line 581

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



700
701
702
703
704
 
# File 'lib/extendmatrix.rb', line 700

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



405
406
407
408
409
410
411
412
 
# File 'lib/extendmatrix.rb', line 405

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



532
533
534
 
# File 'lib/extendmatrix.rb', line 532

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

#sscp ⇒ Object

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



536
537
538
 
# File 'lib/extendmatrix.rb', line 536

def sscp
  transpose*self
end

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

Returns a string for nice printing matrix



458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
 
# File 'lib/extendmatrix.rb', line 458

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 



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

alias :to_s_old :to_s

#total_sum ⇒ Object

Calculate sum of cells



712
713
714
 
# File 'lib/extendmatrix.rb', line 712

def total_sum
  row_sum.inject(&:+)
end

#U ⇒ Object

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



761
762
763
 
# File 'lib/extendmatrix.rb', line 761

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

#wraplate(ijwrap = "") ⇒ Object 



414
415
416
417
418
419
420
421
422
423
424
 
# File 'lib/extendmatrix.rb', line 414

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