Module: NuLin

Defined in:
lib/nulin.rb,
lib/nulin/qr.rb,
lib/nulin/det.rb,
lib/nulin/lls.rb,
lib/nulin/svd.rb,
lib/nulin/gemm.rb,
lib/nulin/cholesky.rb,
lib/nulin/eigensystem.rb,
ext/nulin_native.c

Defined Under Namespace

Modules: Native Classes: Cholesky, ComplexEigenDecomposition, DimensionError, EigenDecomposition, LLS, LLS_QR, LLS_SVD, LinalgError, QR, RealEigenDecomposition, SVD

Constant Summary collapse

TYPECODES = 
[NArray::SFLOAT, NArray::DFLOAT, NArray::SCOMPLEX, NArray::DCOMPLEX]
LAPACK_PREFIX = 
{
  NArray::SFLOAT => "s", NArray::DFLOAT => "d",
  NArray::SCOMPLEX => "c", NArray::DCOMPLEX => "z",
}

Class Method Summary collapse

Class Method Details

.cholesky(matrix, options = {}) ⇒ Object

Note:

This method doesn’t validate the matrix is symmetric/Hermitian. The only upper/lower half of the matrix is used for the computation.

Compute the Cholesky decomposition for a positive definite symmetric/Hermitian ‘matrix`.

Options

  • :type (default :U) - select whether the result matrix is upper(:U) or lower(:L).

If the matrix is not positive definite, NuLin::LinalgError is raised.

Parameters:

  • matrix (NMatrix) —

    target matrix

  • options (Hash) (defaults to: {}) —

    option hash



16
17
18
19
20
21
22
 
# File 'lib/nulin/cholesky.rb', line 16

def cholesky(matrix, options={})
  unless matrix.square?
    raise DimensionError, "Cholesky decomposition is computable for a square matrix"
  end

  return Cholesky.new(matrix, options).result
end

.det(matrix) ⇒ Object

Compute the determinant of the ‘matrix` and return it.

You get an exception NuLin::DimensionError if the matrix is not square The determinant of given matrix is computed using QR factorization.

Parameters:

  • matrix (NMatrix) —

    target matrix



9
10
11
12
13
14
15
16
 
# File 'lib/nulin/det.rb', line 9

def det(matrix)
  unless matrix.square?
    raise DimensionError, "Determinant is computable only on a square matrix"
  end
  n, = matrix.shape
  r = qr(matrix).R
  (0 ... n).inject(1.0){|u, i| u*r[i,i] }
end

.eigensystem(matrix, options = {}) ⇒ Object

Return the eigensystem(a pair of eigenvalues and left/right eigenvectors) of ‘matrix` as an instance of (a subclass of) EigenDecomposition

You can call this with some ‘options` as follows

* :use_complex (default true) Return complex NArray objects
  If you enable this, the return object (EigenDecomposition)
  holds complex NArray objects.
  If you disable this, it is assumed that all eigenvalues are
  real and the return object holds the real NArray objects.
  If you disable this but there is a complex eigenvalue,
  the exception ArgumentError is raised.
  This options is ignored if `matrix` is complex.
* :use_right (default true) - Compute right eigenvectors
* :use_left (default true) - Compute left eigenvectors

Parameters:

  • matrix (NMatrix) —

    target matrix

  • options (Hash) (defaults to: {}) —

    options



21
22
23
24
25
26
27
28
29
30
31
32
33
34
 
# File 'lib/nulin/eigensystem.rb', line 21

def eigensystem(matrix, options={})
  unless matrix.square?
    raise DimensionError, "The eigensystem is computable only for a square matrix"
  end
  
  case
  when matrix.real?
    RealEigenDecomposition.new(matrix, options)
  when matrix.complex?
    ComplexEigenDecomposition.new(matrix, options)
  else
    raise ArgumentError, "The eigensystem is comptable for a complex/float matrix"
  end
end

.linear_least_square(a, b, options = {}) ⇒ Object

Solve the linear least square problem of min_x |a*x - b|.

The solution is returned as an instance of NuLin::LLS.

When n >= m and rank(a) = n, the problem has a unique solution.

When n < m and rank(a) = n, the equation a*x - b = 0 has infinitely many solutions. In this case, we consider the problem of finding the (unique) solution which minimizes |x|.

If rank(a)=min(m, n), we can find the solution of the two problems using QR or LQ factorization.

In the general case, if rank(a) < min(m, n), the above two problems don’t have a unique solution. Therefore we consider the problem of finding x which minimize |a*x-b| and |x|^2. We can solve the problem using complete orthogonal factorization or SVD.

options

  • :algorithm (default :qr) - select algorithm :qr, :svd

  • :rcond (default -1.0) - The threshold to determine effective rank of matrix a. This parameter is used only for svd algorithm. The number of singular values greater than rcond * (the largest singular value) is regarded as the effective rank.

Example # Using QR/LQ algorithm x = NuLin.linear_least_square(a, b).solution

# Using SVD x = NuLin.linear_least_square(a, b, algorithm: :svd, rcond: 0.0001).solution

Parameters:

  • a (NMatrix) —

    a matrix whose shape is [n, m]

  • b (NVector) —

    a vector whose shape is [m]

  • options (Hash) (defaults to: {}) —

    computation options



40
41
42
43
44
45
46
47
 
# File 'lib/nulin/lls.rb', line 40

def linear_least_square(a, b, options={})
  case options.fetch(:algorithm, :qr)
  when :qr
    NuLin::LLS_QR.new(a, b, options)
  when :svd
    NuLin::LLS_SVD.new(a, b, options)
  end
end

.lls(*args) ⇒ Object 



49
 
# File 'lib/nulin/lls.rb', line 49

def lls(*args); linear_least_square(*args); end

.mm_transpose_character(t) ⇒ Object 



38
39
40
41
42
43
44
45
46
47
48
49
 
# File 'lib/nulin/gemm.rb', line 38

def mm_transpose_character(t)
  case t
  when nil, false
    "N"
  when :T, :transpose
    "T"
  when :C, :H, :adjoint
    "C"
  else
    raise ArgumentError, "Unkdnown transposing spec: #{t.inspect}"
  end
end

.multiply_matrix_and_add(alpha, a, b, beta, c, options = {}) ⇒ Object Also known as: mm

Return alpha*a*b + beta*c



5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
 
# File 'lib/nulin/gemm.rb', line 5

def multiply_matrix_and_add(alpha, a, b, beta, c, options={})
  c = c.dup if options.fetch(:overwrite_c, false)
  trans_a = options.fetch(:trans_a, nil)
  trans_b = options.fetch(:trans_b, nil)

  if !trans_a
    k, n = a.shape
  else
    n, k = a.shape
  end
  lda = a.shape[0]

  if !trans_b
    m, kk = b.shape
  else
    kk, m = b.shape
  end
  ldb = b.shape[0]

  mm, nn = c.shape
  ldc = mm
  
  unless k == kk && m == mm && n == nn
    raise DimensionError, "Invalid matrix shape"
  end
  
  NuLin::Native.call(a.typecode, "gemm",
                     mm_transpose_character(trans_b), mm_transpose_character(trans_a),
                     m, n, k, alpha, b, ldb, a, lda, beta, c, ldc)

  c
end

.qr(matrix, options = {}) ⇒ Object

Compute the QR Decomposition of ‘matrix` and return the result as NuLin::QR object.

You can call this with ‘options`, but now the argument is ignored.

Parameters:

  • matrix (NMatrix) —

    The matrix

  • options (Hash) (defaults to: {}) —

    computation options, not used now.



12
13
14
 
# File 'lib/nulin/qr.rb', line 12

def qr(matrix, options={})
  NuLin::QR.new(matrix, options)
end

.svd(matrix, options = {}) ⇒ Object

Compute the Singular Value Decomposition of ‘matrix` and return the result as an instance of SVD

You can call this with some ‘options` as follows

* :use_U (default true) Compute the matrix U
* :use_V (default true) Compute the matrix V and V^t
* :full_matrix (default true)

From the theory of SVD, Any matrix is decomposable into U * sigma * transpose(V) where

* U, V: orthogonal(real) or unitary(complex) matrices
* sigma: diagonal matrix

Parameters:

  • matrix (NMatrix) —

    target matrix

  • options (Hash) (defaults to: {}) —

    computation options



19
20
21
 
# File 'lib/nulin/svd.rb', line 19

def svd(matrix, options={})
  SVD.new(matrix, options)
end

.to_complex_typecode(typecode) ⇒ Object 



27
28
29
30
 
# File 'lib/nulin.rb', line 27

def to_complex_typecode(typecode)
  return NArray::SCOMPLEX if typecode == NArray::SFLOAT
  return NArray::DCOMPLEX 
end

.to_real_typecode(typecode) ⇒ Object 



21
22
23
24
25
 
# File 'lib/nulin.rb', line 21

def to_real_typecode(typecode)
  return NArray::DFLOAT if typecode == NArray::DCOMPLEX
  return NArray::SFLOAT if typecode == NArray::SCOMPLEX
  return typecode
end