Class: Stick::Matrix

Inherits:
Object
  • Object
show all
Includes:
Enumerable, ExceptionForMatrix, Exceptions
Defined in:
lib/stick/matrix/core.rb,
lib/stick/matrix.rb,
lib/stick/matrix/lu.rb,
lib/stick/matrix/givens.rb,
lib/stick/matrix/jacobi.rb,
lib/stick/matrix/exception.rb,
lib/stick/matrix/hessenberg.rb,
lib/stick/matrix/householder.rb

Overview

The Matrix class represents a mathematical matrix, and provides methods for creating special-case matrices (zero, identity, diagonal, singular, vector), operating on them arithmetically and algebraically, and determining their mathematical properties (trace, rank, inverse, determinant).

Note that although matrices should theoretically be rectangular, this is not enforced by the class.

Also note that the determinant of integer matrices may be incorrectly calculated unless you also require 'mathn'. This may be fixed in the future.

Method Catalogue

To create a matrix:

  • Matrix[*rows]

  • Matrix.[](*rows)

  • Matrix.rows(rows, copy = true)

  • Matrix.columns(columns)

  • Matrix.diagonal(*values)

  • Matrix.scalar(n, value)

  • Matrix.scalar(n, value)

  • Matrix.identity(n)

  • Matrix.unit(n)

  • Matrix.I(n)

  • Matrix.zero(n)

  • Matrix.row_vector(row)

  • Matrix.column_vector(column)

To access Matrix elements/columns/rows/submatrices/properties:

  • [](i, j)

  • #row_size

  • #column_size

  • #row(i)

  • #column(j)

  • #collect

  • #map

  • #minor(*param)

Properties of a matrix:

  • #regular?

  • #singular?

  • #square?

Matrix arithmetic:

  • *(m)

  • +(m)

  • -(m)

  • #/(m)

  • #inverse

  • #inv

  • **

Matrix functions:

  • #determinant

  • #det

  • #rank

  • #trace

  • #tr

  • #transpose

  • #t

Conversion to other data types:

  • #coerce(other)

  • #row_vectors

  • #column_vectors

  • #to_a

String representations:

  • #to_s

  • #inspect

Defined Under Namespace

Modules: Exceptions, Givens, Hessenberg, Householder, Jacobi, LU, MMatrix Classes: Scalar

Instance Attribute Summary collapse

Class Method Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(init_method, *argv) ⇒ Matrix

This method is used by the other methods that create matrices, and is of no use to general users.



250
251
252
253
254
255
 
# File 'lib/stick/matrix/core.rb', line 250

def initialize(*argv)
  return initialize_old(*argv) if argv[0].is_a?(Symbol)
  n, m, val = argv; val = 0 if not val
  f = (block_given?)? lambda {|i,j| yield(i, j)} : lambda {|i,j| val}
  init_rows((0...n).collect {|i| (0...m).collect {|j| f.call(i,j)}}, true)
end

Instance Attribute Details

#rows ⇒ Object (readonly)

Returns the value of attribute rows.



217
218
219
 
# File 'lib/stick/matrix.rb', line 217

def rows
  @rows
end

#wrap ⇒ Object

Returns the value of attribute wrap.



217
218
219
 
# File 'lib/stick/matrix.rb', line 217

def wrap
  @wrap
end

Class Method Details

.[](*rows) ⇒ Object

Creates a matrix where each argument is a row.

Matrix[ [25, 93], [-1, 66] ]
   =>  25 93
       -1 66



124
125
126
 
# File 'lib/stick/matrix/core.rb', line 124

def Matrix.[](*rows)
  new(:init_rows, rows, false)
end

.column_vector(column) ⇒ Object

Creates a single-column matrix where the values of that column are as given in column.

Matrix.column_vector([4,5,6])
  => 4
     5
     6



235
236
237
238
239
240
241
242
243
244
 
# File 'lib/stick/matrix/core.rb', line 235

def Matrix.column_vector(column)
  case column
  when Vector
    Matrix.columns([column.to_a])
  when Array
    Matrix.columns([column])
  else
    Matrix.columns([[column]])
  end
end

.columns(columns) ⇒ Object

Creates a matrix using columns as an array of column vectors.

Matrix.columns([[25, 93], [-1, 66]])
   =>  25 -1
       93 66



146
147
148
149
150
151
152
153
154
155
 
# File 'lib/stick/matrix/core.rb', line 146

def Matrix.columns(columns)
  rows = (0 .. columns[0].size - 1).collect {
    |i|
    (0 .. columns.size - 1).collect {
      |j|
      columns[j][i]
    }
  }
  Matrix.rows(rows, false)
end

.diag(*args) ⇒ Object

Creates a matrix with the given matrices as diagonal blocks



332
333
334
335
336
337
338
339
340
341
342
343
344
 
# File 'lib/stick/matrix.rb', line 332

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) 


358
359
360
361
362
363
364
365
 
# File 'lib/stick/matrix.rb', line 358

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

.diagonal(*values) ⇒ Object

Creates a matrix where the diagonal elements are composed of values.

Matrix.diagonal(9, 5, -3)
  =>  9  0  0
      0  5  0
      0  0 -3



164
165
166
167
168
169
170
171
172
173
 
# File 'lib/stick/matrix/core.rb', line 164

def Matrix.diagonal(*values)
  size = values.size
  rows = (0 .. size  - 1).collect {
    |j|
    row = Array.new(size).fill(0, 0, size)
    row[j] = values[j]
    row
  }
  rows(rows, false)
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) 


349
350
351
352
353
 
# File 'lib/stick/matrix.rb', line 349

def equal_in_delta?(m0, m1, delta = 1.0e-10)
  delta = delta.abs
  mapcar(m0, m1){|x, y| return false if (x < y - delta or x > y + delta)  }
  true
end

.identity(n) ⇒ Object Also known as: unit, I

Creates an n by n identity matrix.

Matrix.identity(2)
  => 1 0
     0 1



192
193
194
 
# File 'lib/stick/matrix/core.rb', line 192

def Matrix.identity(n)
  Matrix.scalar(n, 1)
end

.row_vector(row) ⇒ Object

Creates a single-row matrix where the values of that row are as given in row.

Matrix.row_vector([4,5,6])
  => 4 5 6



216
217
218
219
220
221
222
223
224
225
 
# File 'lib/stick/matrix/core.rb', line 216

def Matrix.row_vector(row)
  case row
  when Vector
    Matrix.rows([row.to_a], false)
  when Array
    Matrix.rows([row.dup], false)
  else
    Matrix.rows([[row]], false)
  end
end

.rows(rows, copy = true) ⇒ Object

Creates a matrix where rows is an array of arrays, each of which is a row to the matrix. If the optional argument copy is false, use the given arrays as the internal structure of the matrix without copying.

Matrix.rows([[25, 93], [-1, 66]])
   =>  25 93
       -1 66



135
136
137
 
# File 'lib/stick/matrix/core.rb', line 135

def Matrix.rows(rows, copy = true)
  new(:init_rows, rows, copy)
end

.scalar(n, value) ⇒ Object

Creates an n by n diagonal matrix where each diagonal element is value.

Matrix.scalar(2, 5)
  => 5 0
     0 5



182
183
184
 
# File 'lib/stick/matrix/core.rb', line 182

def Matrix.scalar(n, value)
  Matrix.diagonal(*Array.new(n).fill(value, 0, n))
end

.zero(n) ⇒ Object

Creates an n by n zero matrix.

Matrix.zero(2)
  => 0 0
     0 0



206
207
208
 
# File 'lib/stick/matrix/core.rb', line 206

def Matrix.zero(n)
  Matrix.scalar(n, 0)
end

Instance Method Details

#*(m) ⇒ Object

Matrix multiplication.

Matrix[[2,4], [6,8]] * Matrix.identity(2)
  => 2 4
     6 8



462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
 
# File 'lib/stick/matrix/core.rb', line 462

def *(m) # m is matrix or vector or number
  case(m)
  when Numeric
    rows = @rows.collect {
      |row|
      row.collect {
        |e|
        e * m
      }
    }
    return Matrix.rows(rows, false)
  when Vector
    m = Matrix.column_vector(m)
    r = self * m
    return r.column(0)
  when Matrix
    Matrix.Raise ErrDimensionMismatch if column_size != m.row_size

    rows = (0 .. row_size - 1).collect {
      |i|
      (0 .. m.column_size - 1).collect {
        |j|
        vij = 0
        0.upto(column_size - 1) do
          |k|
          vij += self[i, k] * m[k, j]
        end
        vij
      }
    }
    return Matrix.rows(rows, false)
  else
    x, y = m.coerce(self)
    return x * y
  end
end

#**(other) ⇒ Object

Matrix exponentiation. Defined for integer powers only. Equivalent to multiplying the matrix by itself N times.

Matrix[[7,6], [3,9]] ** 2
  => 67 96
     48 99



655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
 
# File 'lib/stick/matrix/core.rb', line 655

def ** (other)
  if other.kind_of?(Integer)
    x = self
    if other <= 0
      x = self.inverse
      return Matrix.identity(self.column_size) if other == 0
      other = -other
    end
    z = x
    n = other  - 1
    while n != 0
      while (div, mod = n.divmod(2)
             mod == 0)
        x = x * x
        n = div
      end
      z *= x
      n -= 1
    end
    z
  elsif other.kind_of?(Float) || defined?(Rational) && other.kind_of?(Rational)
    Matrix.Raise ErrOperationNotDefined, "**"
  else
    Matrix.Raise ErrOperationNotDefined, "**"
  end
end

#+(m) ⇒ Object

Matrix addition.

Matrix.scalar(2,5) + Matrix[[1,0], [-4,7]]
  =>  6  0
     -4 12



505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
 
# File 'lib/stick/matrix/core.rb', line 505

def +(m)
  case m
  when Numeric
    Matrix.Raise ErrOperationNotDefined, "+"
  when Vector
    m = Matrix.column_vector(m)
  when Matrix
  else
    x, y = m.coerce(self)
    return x + y
  end

  Matrix.Raise ErrDimensionMismatch unless row_size == m.row_size and column_size == m.column_size

  rows = (0 .. row_size - 1).collect {
    |i|
    (0 .. column_size - 1).collect {
      |j|
      self[i, j] + m[i, j]
    }
  }
  Matrix.rows(rows, false)
end

#-(m) ⇒ Object

Matrix subtraction.

Matrix[[1,5], [4,2]] - Matrix[[9,3], [-4,1]]
  => -8  2
      8  1



535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
 
# File 'lib/stick/matrix/core.rb', line 535

def -(m)
  case m
  when Numeric
    Matrix.Raise ErrOperationNotDefined, "-"
  when Vector
    m = Matrix.column_vector(m)
  when Matrix
  else
    x, y = m.coerce(self)
    return x - y
  end

  Matrix.Raise ErrDimensionMismatch unless row_size == m.row_size and column_size == m.column_size

  rows = (0 .. row_size - 1).collect {
    |i|
    (0 .. column_size - 1).collect {
      |j|
      self[i, j] - m[i, j]
    }
  }
  Matrix.rows(rows, false)
end

#==(other) ⇒ Object Also known as: eql?

Returns true if and only if the two matrices contain equal elements.



411
412
413
414
415
 
# File 'lib/stick/matrix/core.rb', line 411

def ==(other)
  return false unless Matrix === other

  other.compare_by_row_vectors(@rows)
end

#[](i, j) ⇒ Object Also known as: element, component

Returns element (i,j) of the matrix. That is: row i, column j.



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

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 Also known as: set_element, set_component

Set the values of a matrix m = Matrix.new(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]]



283
284
285
286
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
 
# File 'lib/stick/matrix.rb', line 283

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



79
80
81
82
 
# File 'lib/stick/matrix/householder.rb', line 79

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

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

Classical Jacobi 8.4.3 Golub & van Loan



66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
 
# File 'lib/stick/matrix/jacobi.rb', line 66

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.



85
86
87
 
# File 'lib/stick/matrix/jacobi.rb', line 85

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.



100
101
102
 
# File 'lib/stick/matrix/jacobi.rb', line 100

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

#clone ⇒ Object

Returns a clone of the matrix, so that the contents of each do not reference identical objects.



318
319
320
 
# File 'lib/stick/matrix.rb', line 318

def clone
  super
end

#coerce(other) ⇒ Object

FIXME: describe #coerce.



910
911
912
913
914
915
916
917
 
# File 'lib/stick/matrix/core.rb', line 910

def coerce(other)
  case other
  when Numeric
    return Scalar.new(other), self
  else
    raise TypeError, "#{self.class} can't be coerced into #{other.class}"
  end
end

#collect ⇒ Object Also known as: map

Returns a matrix that is the result of iteration of the given block over all elements of the matrix.

Matrix[ [1,2], [3,4] ].collect { |e| e**2 }
  => 1  4
     9 16



338
339
340
341
 
# File 'lib/stick/matrix/core.rb', line 338

def collect # :yield: e
  rows = @rows.collect{|row| row.collect{|e| yield e}}
  Matrix.rows(rows, false)
end

#cols_len ⇒ Object

Returns a list with the maximum lengths



427
428
429
 
# File 'lib/stick/matrix.rb', line 427

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

#column(j) ⇒ Object

Returns column vector number j of the matrix as a Vector (starting at 0 like an array). When a block is given, the elements of that vector are iterated.



316
317
318
319
320
321
322
323
324
325
326
327
328
329
 
# File 'lib/stick/matrix/core.rb', line 316

def column(j) # :yield: e
  if block_given?
    0.upto(row_size - 1) do
      |i|
      yield @rows[i][j]
    end
  else
    col = (0 .. row_size - 1).collect {
      |i|
      @rows[i][j]
    }
    Vector.elements(col, false)
  end
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



520
521
522
523
524
525
526
 
# File 'lib/stick/matrix.rb', line 520

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 colomn/s of matrix as a Matrix



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

def column2matrix(c)
  a = self.send(:column, c).to_a
  if c.is_a?(Range) and c.entries.size > 1
    return Matrix[*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.new(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



537
538
539
540
541
542
543
 
# File 'lib/stick/matrix.rb', line 537

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.



511
512
513
514
 
# File 'lib/stick/matrix.rb', line 511

def column_collect(j, &block)
  f = MMatrix.default_block(block)
  (0...row_size).collect {|r| f.call(self[r, j])}
end

#column_size ⇒ Object

Returns the number of columns. Note that it is possible to construct a matrix with uneven columns (e.g. Matrix[ [1,2,3], [4,5] ]), but this is mathematically unsound. This method uses the first row to determine the result.



293
294
295
 
# File 'lib/stick/matrix/core.rb', line 293

def column_size
  @rows[0].size
end

#column_vectors ⇒ Object

Returns an array of the column vectors of the matrix. See Vector.



933
934
935
936
937
938
939
 
# File 'lib/stick/matrix/core.rb', line 933

def column_vectors
  columns = (0 .. column_size - 1).collect {
    |i|
    column(i)
  }
  columns
end

#compare_by_row_vectors(rows) ⇒ Object

Not really intended for general consumption.



421
422
423
424
425
426
427
428
429
 
# File 'lib/stick/matrix/core.rb', line 421

def compare_by_row_vectors(rows)
  return false unless @rows.size == rows.size

  0.upto(@rows.size - 1) do
    |i|
    return false unless @rows[i] == rows[i]
  end
  true
end

#determinant ⇒ Object Also known as: det

Returns the determinant of the matrix. If the matrix is not square, the result is 0. This method’s algorism is Gaussian elimination method and using Numeric#quo(). Beware that using Float values, with their usual lack of precision, can affect the value returned by this method. Use Rational values or Matrix#det_e instead if this is important to you.

Matrix[[7,6], [3,9]].determinant
  => 63.0



696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
 
# File 'lib/stick/matrix/core.rb', line 696

def determinant
  return 0 unless square?

  size = row_size - 1
  a = to_a

  det = 1
  k = 0
  begin
    if (akk = a[k][k]) == 0
      i = k
      begin
        return 0 if (i += 1) > size
      end while a[i][k] == 0
      a[i], a[k] = a[k], a[i]
      akk = a[k][k]
      det *= -1
    end
    (k + 1).upto(size) do
      |i|
      q = a[i][k].quo(akk)
      (k + 1).upto(size) do
        |j|
        a[i][j] -= a[k][j] * q
      end
    end
    det *= akk
  end while (k += 1) <= size
  det
end

#determinant_e ⇒ Object Also known as: det_e

Returns the determinant of the matrix. If the matrix is not square, the result is 0. This method’s algorism is Gaussian elimination method. This method uses Euclidean algorism. If all elements are integer, really exact value. But, if an element is a float, can’t return exact value.

Matrix[[7,6], [3,9]].determinant
  => 63



738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
 
# File 'lib/stick/matrix/core.rb', line 738

def determinant_e
  return 0 unless square?

  size = row_size - 1
  a = to_a

  det = 1
  k = 0
  begin
    if a[k][k].zero?
      i = k
      begin
        return 0 if (i += 1) > size
      end while a[i][k].zero?
      a[i], a[k] = a[k], a[i]
      det *= -1
    end
    (k + 1).upto(size) do |i|
      q = a[i][k].quo(a[k][k])
      k.upto(size) do |j|
        a[i][j] -= a[k][j] * q
      end
      unless a[i][k].zero?
        a[i], a[k] = a[k], a[i]
        det *= -1
        redo
      end
    end
    det *= a[k][k]
  end while (k += 1) <= size
  det
end

#each ⇒ Object

Iterate the elements of a matrix



453
454
455
456
 
# File 'lib/stick/matrix.rb', line 453

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

#eigenvaluesJacobi ⇒ Object

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



92
93
94
95
 
# File 'lib/stick/matrix/jacobi.rb', line 92

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

#elements_to_f ⇒ Object 



948
949
950
 
# File 'lib/stick/matrix/core.rb', line 948

def elements_to_f
  collect{|e| e.to_f}
end

#elements_to_i ⇒ Object 



952
953
954
 
# File 'lib/stick/matrix/core.rb', line 952

def elements_to_i
  collect{|e| e.to_i}
end

#elements_to_r ⇒ Object 



956
957
958
 
# File 'lib/stick/matrix/core.rb', line 956

def elements_to_r
  collect{|e| e.to_r}
end

#empty? ⇒ Boolean

Tests if the matrix is empty or not

Returns:

  • (Boolean) 


570
571
572
 
# File 'lib/stick/matrix.rb', line 570

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



53
54
55
 
# File 'lib/stick/matrix/givens.rb', line 53

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

#givensR ⇒ Object

Returns the upper triunghiular matrix R of a Givens QR factorization



45
46
47
 
# File 'lib/stick/matrix/givens.rb', line 45

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

#gram_schmidt ⇒ Object

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



603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
 
# File 'lib/stick/matrix.rb', line 603

def gram_schmidt
  a = clone
  n = column_size
  m = row_size
  q = Matrix.new(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



624
625
626
 
# File 'lib/stick/matrix.rb', line 624

def gram_schmidtQ
  gram_schmidt[0]
end

#gram_schmidtR ⇒ Object

Returns the R_1 upper triangular matrix of Modified Gram Schmidt algorithm



631
632
633
 
# File 'lib/stick/matrix.rb', line 631

def gram_schmidtR
  gram_schmidt[1]
end

#hash ⇒ Object

Returns a hash-code for the matrix.



442
443
444
445
446
447
448
449
450
 
# File 'lib/stick/matrix/core.rb', line 442

def hash
  value = 0
  for row in @rows
    for e in row
      value ^= e.hash
    end
  end
  return value
end

#hessenberg_form_H ⇒ Object

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



40
41
42
 
# File 'lib/stick/matrix/hessenberg.rb', line 40

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)



28
29
30
 
# File 'lib/stick/matrix/hessenberg.rb', line 28

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

#hessenbergR ⇒ Object

Returns the upper triunghiular matrix R of a Hessenberg QR factorization



34
35
36
 
# File 'lib/stick/matrix/hessenberg.rb', line 34

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,



87
88
89
90
91
92
 
# File 'lib/stick/matrix/householder.rb', line 87

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



97
98
99
100
101
102
 
# File 'lib/stick/matrix/householder.rb', line 97

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

#ids ⇒ Object 



240
 
# File 'lib/stick/matrix.rb', line 240

alias :ids :[]

#initialize_copy(orig) ⇒ Object 



322
323
324
325
 
# File 'lib/stick/matrix.rb', line 322

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

#initialize_old(init_method, *argv) ⇒ Object

For invoking a method in Ruby1.8 is working ‘send’ and in Ruby1.9 is working ‘funcall’



232
233
234
235
236
237
238
 
# File 'lib/stick/matrix.rb', line 232

def initialize_old(init_method, *argv)
  if RUBY_VERSION < "1.9.0"
    self.send(init_method, *argv) # in Ruby1.8
  else
    self.funcall(init_method, *argv) # in Ruby1.9
  end
end

#inspect ⇒ Object

Overrides Object#inspect



977
978
979
 
# File 'lib/stick/matrix/core.rb', line 977

def inspect
  "Matrix"+@rows.inspect
end

#inverse ⇒ Object Also known as: inv

Returns the inverse of the matrix.

Matrix[[1, 2], [2, 1]].inverse
  => -1  1
      0 -1



590
591
592
593
 
# File 'lib/stick/matrix/core.rb', line 590

def inverse
  Matrix.Raise ErrDimensionMismatch unless square?
  Matrix.I(row_size).inverse_from(self)
end

#inverse_from(src) ⇒ Object

Not for public consumption?



599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
 
# File 'lib/stick/matrix/core.rb', line 599

def inverse_from(src)
  size = row_size - 1
  a = src.to_a

  for k in 0..size
    i = k
    akk = a[k][k].abs
    for j in (k+1)..size
      v = a[j][k].abs
      if v > akk
        i = j
        akk = v
      end
    end
    Matrix.Raise ErrNotRegular if akk == 0
    if i != k
      a[i], a[k] = a[k], a[i]
      @rows[i], @rows[k] = @rows[k], @rows[i]
    end
    akk = a[k][k]

    for i in 0 .. size
      next if i == k
      q = a[i][k].quo(akk)
      a[i][k] = 0

      (k + 1).upto(size) do
        |j|
        a[i][j] -= a[k][j] * q
      end
      0.upto(size) do
        |j|
        @rows[i][j] -= @rows[k][j] * q
      end
    end

    (k + 1).upto(size) do
      |j|
      a[k][j] = a[k][j].quo(akk)
    end
    0.upto(size) do
      |j|
      @rows[k][j] = @rows[k][j].quo(akk)
    end
  end
  self
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



54
55
56
 
# File 'lib/stick/matrix/lu.rb', line 54

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

#max_len_column(j) ⇒ Object

Returns the maximum length of column elements



420
421
422
 
# File 'lib/stick/matrix.rb', line 420

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

#minor(*param) ⇒ Object

Returns a section of the matrix. The parameters are either:

  • start_row, nrows, start_col, ncols; OR

  • col_range, row_range

Matrix.diagonal(9, 5, -3).minor(0..1, 0..2)
  => 9 0 0
     0 5 0



353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
 
# File 'lib/stick/matrix/core.rb', line 353

def minor(*param)
  case param.size
  when 2
    from_row = param[0].first
    size_row = param[0].end - from_row
    size_row += 1 unless param[0].exclude_end?
    from_col = param[1].first
    size_col = param[1].end - from_col
    size_col += 1 unless param[1].exclude_end?
  when 4
    from_row = param[0]
    size_row = param[1]
    from_col = param[2]
    size_col = param[3]
  else
    Matrix.Raise ArgumentError, param.inspect
  end

  rows = @rows[from_row, size_row].collect{
    |row|
    row[from_col, size_col]
  }
  Matrix.rows(rows, false)
end

#norm(p = 2) ⇒ Object 



558
559
560
 
# File 'lib/stick/matrix.rb', line 558

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

#norm_frobenius ⇒ Object Also known as: normF 



562
563
564
 
# File 'lib/stick/matrix.rb', line 562

def norm_frobenius
  norm
end

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

Returns the matrix divided by a scalar



371
372
373
 
# File 'lib/stick/matrix.rb', line 371

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

#rank ⇒ Object

Returns the rank of the matrix. Beware that using Float values, probably return faild value. Use Rational values or Matrix#rank_e for getting exact result.

Matrix[[7,6], [3,9]].rank
  => 2



780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
 
# File 'lib/stick/matrix/core.rb', line 780

def rank
  if column_size > row_size
    a = transpose.to_a
    a_column_size = row_size
    a_row_size = column_size
  else
    a = to_a
    a_column_size = column_size
    a_row_size = row_size
  end
  rank = 0
  k = 0
  begin
    if (akk = a[k][k]) == 0
      i = k
      exists = true
      begin
        if (i += 1) > a_column_size - 1
          exists = false
          break
        end
      end while a[i][k] == 0
      if exists
        a[i], a[k] = a[k], a[i]
        akk = a[k][k]
      else
        i = k
        exists = true
        begin
          if (i += 1) > a_row_size - 1
            exists = false
            break
          end
        end while a[k][i] == 0
        if exists
          k.upto(a_column_size - 1) do
            |j|
            a[j][k], a[j][i] = a[j][i], a[j][k]
          end
          akk = a[k][k]
        else
          next
        end
      end
    end
    (k + 1).upto(a_row_size - 1) do
      |i|
      q = a[i][k].quo(akk)
      (k + 1).upto(a_column_size - 1) do
        |j|
        a[i][j] -= a[k][j] * q
      end
    end
    rank += 1
  end while (k += 1) <= a_column_size - 1
  return rank
end

#rank_e ⇒ Object

Returns the rank of the matrix. This method uses Euclidean algorism. If all elements are integer, really exact value. But, if an element is a float, can’t return exact value.

Matrix[[7,6], [3,9]].rank
  => 2



846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
 
# File 'lib/stick/matrix/core.rb', line 846

def rank_e
  a = to_a
  a_column_size = column_size
  a_row_size = row_size
  pi = 0
  (0 ... a_column_size).each do |j|
    if i = (pi ... a_row_size).find{|i0| !a[i0][j].zero?}
      if i != pi
        a[pi], a[i] = a[i], a[pi]
      end
      (pi + 1 ... a_row_size).each do |k|
        q = a[k][j].quo(a[pi][j])
        (pi ... a_column_size).each do |j0|
          a[k][j0] -= q * a[pi][j0]
        end
        if k > pi && !a[k][j].zero?
          a[k], a[pi] = a[pi], a[k]
          redo
        end
      end
      pi += 1
    end
  end
  pi
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



47
48
49
50
51
52
53
54
55
56
57
58
59
 
# File 'lib/stick/matrix/hessenberg.rb', line 47

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

#regular? ⇒ Boolean

Returns true if this is a regular matrix.

Returns:

  • (Boolean) 


385
386
387
 
# File 'lib/stick/matrix/core.rb', line 385

def regular?
  square? and rank == column_size
end

#row(i) ⇒ Object

Returns row vector number i of the matrix as a Vector (starting at 0 like an array). When a block is given, the elements of that vector are iterated.



301
302
303
304
305
306
307
308
309
 
# File 'lib/stick/matrix/core.rb', line 301

def row(i) # :yield: e
  if block_given?
    for e in @rows[i]
      yield e
    end
  else
    Vector.elements(@rows[i])
  end
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



498
499
500
501
502
503
504
 
# File 'lib/stick/matrix.rb', line 498

def row!(i)
  if block_given?
    @rows[i].collect! {|e| yield e }
  else
    Vector.elements(@rows[i], false)
  end
end

#row2matrix(r) ⇒ Object

Returns the row/s of matrix as a Matrix



577
578
579
580
581
582
583
584
 
# File 'lib/stick/matrix.rb', line 577

def row2matrix(r)
  a = self.send(:row, r).to_a
  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



553
554
555
556
 
# File 'lib/stick/matrix.rb', line 553

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.



489
490
491
492
 
# File 'lib/stick/matrix.rb', line 489

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

#row_size ⇒ Object

Returns the number of rows.



283
284
285
 
# File 'lib/stick/matrix/core.rb', line 283

def row_size
  @rows.size
end

#row_vectors ⇒ Object

Returns an array of the row vectors of the matrix. See Vector.



922
923
924
925
926
927
928
 
# File 'lib/stick/matrix/core.rb', line 922

def row_vectors
  rows = (0 .. row_size - 1).collect {
    |i|
    row(i)
  }
  rows
end

#set(m) ⇒ Object

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



383
384
385
386
387
388
389
390
 
# File 'lib/stick/matrix.rb', line 383

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

#singular? ⇒ Boolean

Returns true is this is a singular (i.e. non-regular) matrix.

Returns:

  • (Boolean) 


392
393
394
 
# File 'lib/stick/matrix/core.rb', line 392

def singular?
  not regular?
end

#square? ⇒ Boolean

Returns true is this is a square matrix. See note in column_size about this being unreliable, though.

Returns:

  • (Boolean) 


400
401
402
 
# File 'lib/stick/matrix/core.rb', line 400

def square?
  column_size == row_size
end

#to_a ⇒ Object

Returns an array of arrays that describe the rows of the matrix.



944
945
946
 
# File 'lib/stick/matrix/core.rb', line 944

def to_a
  @rows.collect{|row| row.collect{|e| e}}
end

#to_s ⇒ Object

Overrides Object#to_s



434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
 
# File 'lib/stick/matrix.rb', line 434

def to_s(mode = :pretty, len_col = 3)
  return super if empty?
  if mode == :pretty
    clen = cols_len
    to_a.collect {|r| mapcar(r, clen) {|x, l| 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

#trace ⇒ Object Also known as: tr

Returns the trace (sum of diagonal elements) of the matrix.

Matrix[[7,6], [3,9]].trace
  => 16



878
879
880
881
882
883
884
885
 
# File 'lib/stick/matrix/core.rb', line 878

def trace
  tr = 0
  0.upto(column_size - 1) do
    |i|
    tr += @rows[i][i]
  end
  tr
end

#transpose ⇒ Object Also known as: t

Returns the transpose of the matrix.

Matrix[[1,2], [3,4], [5,6]]
  => 1 2
     3 4
     5 6
Matrix[[1,2], [3,4], [5,6]].transpose
  => 1 3 5
     2 4 6



898
899
900
 
# File 'lib/stick/matrix/core.rb', line 898

def transpose
  Matrix.columns(@rows)
end

#U ⇒ Object

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



47
48
49
 
# File 'lib/stick/matrix/lu.rb', line 47

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

#wraplate(ijwrap = "") ⇒ Object 



392
393
394
395
396
397
398
399
400
401
402
 
# File 'lib/stick/matrix.rb', line 392

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