Module: Numo::Linalg

Defined in:
lib/numo/linalg/loader.rb,
lib/numo/linalg/version.rb,
lib/numo/linalg/function.rb,
lib/numo/linalg/autoloader.rb,
ext/numo/linalg/blas/blas.c,
ext/numo/linalg/lapack/lapack.c

Defined Under Namespace

Modules: Autoloader, Blas, Lapack, Loader

Constant Summary collapse

VERSION = 
"0.1.7"
BLAS_CHAR = 
{
 SFloat => "s",
 DFloat => "d",
 SComplex => "c",
 DComplex => "z",
}

Class Method Summary collapse

Class Method Details

.blas_char(*args) ⇒ Object 



71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
 
# File 'lib/numo/linalg/function.rb', line 71

def blas_char(*args)
  t = Float
  args.each do |a|
    k =
      case a
      when NArray
        a.class
      when Array
        NArray.array_type(a)
      end
    if k && k < NArray
      t = k::UPCAST[t] || t::UPCAST[k]
    end
  end
  BLAS_CHAR[t] || raise(TypeError,"invalid data type for BLAS/LAPACK")
end

.cho_fact(a, uplo: 'U') ⇒ Numo::NArray

Computes the Cholesky factorization of a symmetric/Hermitian positive definite matrix A. The factorization has the form

A = U**H * U,  if UPLO = 'U', or
A = L  * L**H,  if UPLO = 'L',

where U is an upper triangular matrix and L is a lower triangular matrix.

Parameters:

  • a (Numo::NArray) —

    n-by-n symmetric matrix A (>= 2-dimensinal NArray)

  • uplo (String or Symbol) (defaults to: 'U') —

    optional, default=‘U’. Access upper or (‘U’) lower (‘L’) triangle.

Returns:

  • (Numo::NArray) —

    The matrix which has the Cholesky factor in upper or lower triangular part. Remain part consists of random values.



591
592
593
 
# File 'lib/numo/linalg/function.rb', line 591

def cho_fact(a, uplo:'U')
  Lapack.call(:potrf, a, uplo:uplo)[0]
end

.cho_inv(a, uplo: 'U') ⇒ Numo::NArray

Computes the inverse of a symmetric/Hermitian positive definite matrix A using the Cholesky factorization ‘A = U**T*U` or `A = L*L**T` computed by Linalg.cho_fact.

Parameters:

  • a (Numo::NArray) —

    the triangular factor U or L from the Cholesky factorization, as computed by Linalg.cho_fact.

  • uplo (String or Symbol) (defaults to: 'U') —

    optional, default=‘U’. Access upper or (‘U’) lower (‘L’) triangle.

Returns:

  • (Numo::NArray) —

    the upper or lower triangle of the (symmetric) inverse of A.



606
607
608
 
# File 'lib/numo/linalg/function.rb', line 606

def cho_inv(a, uplo:'U')
  Lapack.call(:potri, a, uplo:uplo)[0]
end

.cho_solve(a, b, uplo: 'U') ⇒ Numo::NArray

Solves a system of linear equations

A*X = B

with a symmetric/Hermitian positive definite matrix A using the Cholesky factorization ‘A = U**T*U` or `A = L*L**T` computed by Linalg.cho_fact.

Parameters:

  • a (Numo::NArray) —

    the triangular factor U or L from the Cholesky factorization, as computed by Linalg.cho_fact.

  • b (Numo::NArray) —

    the right hand side matrix B.

  • uplo (String or Symbol) (defaults to: 'U') —

    optional, default=‘U’. Access upper or (‘U’) lower (‘L’) triangle.

Returns:

  • (Numo::NArray) —

    the solution matrix X.



622
623
624
 
# File 'lib/numo/linalg/function.rb', line 622

def cho_solve(a, b, uplo:'U')
  Lapack.call(:potrs, a, b, uplo:uplo)[0]
end

.cholesky(a, uplo: 'U') ⇒ Numo::NArray

Computes the Cholesky factorization of a symmetric/Hermitian positive definite matrix A. The factorization has the form

A = U**H * U,  if UPLO = 'U', or
A = L  * L**H,  if UPLO = 'L',

where U is an upper triangular matrix and L is a lower triangular matrix.

Parameters:

  • a (Numo::NArray) —

    n-by-n symmetric matrix A (>= 2-dimensinal NArray)

  • uplo (String or Symbol) (defaults to: 'U') —

    optional, default=‘U’. Access upper or (‘U’) lower (‘L’) triangle.

Returns:

  • (Numo::NArray) —

    The factor U or L.

Raises:

  • (NArray::ShapeError) 


564
565
566
567
568
569
570
571
572
573
574
575
576
 
# File 'lib/numo/linalg/function.rb', line 564

def cholesky(a, uplo: 'U')
  raise NArray::ShapeError, '2-d array is required' if a.ndim < 2
  raise NArray::ShapeError, 'matrix a is not square matrix' if a.shape[0] != a.shape[1]
  factor = Lapack.call(:potrf, a, uplo: uplo)[0]
  if uplo == 'U'
    factor.triu
  else
    # TODO: Use the tril method if the verision of Numo::NArray
    #       in the runtime dependency list becomes 0.9.1.3 or higher.
    m, = a.shape
    factor * Numo::DFloat.ones(m, m).triu.transpose
  end
end

.cond(a, ord = nil) ⇒ Numo::NArray

Compute the condition number of a matrix using the norm with one of the following order.

|  ord  |  matrix norm           |
| ----- | ---------------------- |
|  nil  | 2-norm using SVD       |
| 'fro' | Frobenius norm         |
| 'inf' | x.abs.sum(axis:-1).max |
|    1  | x.abs.sum(axis:-2).max |
|    2  | 2-norm (max sing_vals) |

Examples:

a = Numo::DFloat[[1, 0, -1], [0, 1, 0], [1, 0, 1]]
=> Numo::DFloat#shape=[3,3]
[[1, 0, -1],
 [0, 1, 0],
 [1, 0, 1]]
LA = Numo::Linalg
LA.cond(a)
=> 1.4142135623730951
LA.cond(a, 'fro')
=> 3.1622776601683795
LA.cond(a, 'inf')
=> 2.0
LA.cond(a, '-inf')
=> 1.0
LA.cond(a, 1)
=> 2.0
LA.cond(a, -1)
=> 1.0
LA.cond(a, 2)
=> 1.4142135623730951
LA.cond(a, -2)
=> 0.7071067811865475
(LA.svdvals(a)).min*(LA.svdvals(LA.inv(a))).min
=> 0.7071067811865475

Parameters:

  • a (Numo::NArray) —

    matrix or vector (>= 1-dimensinal NArray)

  • ord (String or Symbol) (defaults to: nil) —

    Order of the norm.

Returns:

  • (Numo::NArray) —

    cond result



896
897
898
899
900
901
902
903
 
# File 'lib/numo/linalg/function.rb', line 896

def cond(a,ord=nil)
  if ord.nil?
    s = svdvals(a)
    s[false, 0]/s[false, -1]
  else
    norm(a, ord, axis:[-2,-1]) * norm(inv(a), ord, axis:[-2,-1])
  end
end

.det(a) ⇒ Float or Complex or Numo::NArray

Determinant of a matrix

Parameters:

  • a (Numo::NArray) —

    matrix (>= 2-dimensional NArray)

Returns:

  • (Float or Complex or Numo::NArray) 


910
911
912
913
914
915
916
 
# File 'lib/numo/linalg/function.rb', line 910

def det(a)
  lu, piv, = Lapack.call(:getrf, a)
  idx = piv.new_narray.store(piv.class.new(piv.shape[-1]).seq(1))
  m = piv.eq(idx).count_false(axis:-1) % 2
  sign = m * -2 + 1
  lu.diagonal.prod(axis:-1) * sign
end

.dot(a, b) ⇒ Numo::NArray

Dot product.

Parameters:

  • a (Numo::NArray) —

    matrix or vector (>= 1-dimensinal NArray)

  • b (Numo::NArray) —

    matrix or vector (>= 1-dimensinal NArray)

Returns:

  • (Numo::NArray) —

    result of dot product



96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
 
# File 'lib/numo/linalg/function.rb', line 96

def dot(a, b)
  a = NArray.asarray(a)
  b = NArray.asarray(b)
  case a.ndim
  when 1
    case b.ndim
    when 1
      func = blas_char(a, b) =~ /c|z/ ? :dotu : :dot
      Blas.call(func, a, b)
    else
      if b.contiguous?
        trans = 't'
      else
        if b.fortran_contiguous?
          trans = 'n'
          b = b.transpose
        else
          trans = 't'
          b = b.dup
        end
      end
      Blas.call(:gemv, b, a, trans:trans)
    end
  else
    case b.ndim
    when 1
      if a.contiguous?
        trans = 'n'
      else
        if a.fortran_contiguous?
          trans = 't'
          a = a.transpose
        else
          trans = 'n'
          a = a.dup
        end
      end
      Blas.call(:gemv, a, b, trans:trans)
    else
      if a.contiguous?
        transa = 'n'
      else
        if a.fortran_contiguous?
          transa = 't'
          a = a.transpose
        else
          transa = 'n'
          a = a.dup
        end
      end
      if b.contiguous?
        transb = 'n'
      else
        if b.fortran_contiguous?
          transb='t'
          b = b.transpose
        else
          transb='n'
          b = b.dup
        end
      end
      Blas.call(:gemm, a, b, transa:transa, transb:transb)
    end
  end
end

.eig(a, left: false, right: true) ⇒ [w,vl,vr]

Computes the eigenvalues and, optionally, the left and/or right eigenvectors for a square nonsymmetric matrix A.

Parameters:

  • a (Numo::NArray) —

    square nonsymmetric matrix (>= 2-dimensinal NArray)

  • left (Bool) (defaults to: false) —

    (optional) If true, left eigenvectors are computed.

  • right (Bool) (defaults to: true) —

    (optional) If true, right eigenvectors are computed.

Returns:

  • ([w,vl,vr]) —

    • w [Numo::NArray] – The eigenvalues.

    • vl [Numo::NArray] – The left eigenvectors if left is true, otherwise nil.

    • vr [Numo::NArray] – The right eigenvectors if right is true, otherwise nil.



640
641
642
643
644
645
646
647
648
649
650
651
652
 
# File 'lib/numo/linalg/function.rb', line 640

def eig(a, left:false, right:true)
  jobvl, jobvr = left, right
  case blas_char(a)
  when /c|z/
    w, vl, vr, info = Lapack.call(:geev, a, jobvl:jobvl, jobvr:jobvr)
  else
    wr, wi, vl, vr, info = Lapack.call(:geev, a, jobvl:jobvl, jobvr:jobvr)
    w  = wr + wi * Complex::I
    vl = _make_complex_eigvecs(w,vl) if left
    vr = _make_complex_eigvecs(w,vr) if right
  end
  [w,vl,vr] #.compact
end

.eigh(a, b = nil, vals_only: false, vals_range: nil, uplo: 'U', turbo: false) ⇒ [w,v]

Obtains the eigenvalues and, optionally, the eigenvectors by solving an ordinary or generalized eigenvalue problem for a square symmetric / Hermitian matrix.

Parameters:

  • a (Numo::NArray) —

    square symmetric matrix (>= 2-dimensinal NArray)

  • b (Numo::NArray) (defaults to: nil) —

    (optional) square symmetric matrix (>= 2-dimensinal NArray) If nil, identity matrix is assumed.

  • vals_only (Bool) (defaults to: false) —

    (optional) If false, eigenvectors are computed.

  • vals_range (Range) (defaults to: nil) —

    (optional) The range of indices of the eigenvalues (in ascending order) and corresponding eigenvectors to be returned. If nil or 0…N (N is the size of the matrix a), all eigenvalues and eigenvectors are returned.

  • uplo (String or Symbol) (defaults to: 'U') —

    (optional, default=‘U’) Access upper (‘U’) or lower (‘L’) triangle.

  • turbo (Bool) (defaults to: false) —

    (optional) If true, divide and conquer algorithm is used.

Returns:

  • ([w,v]) —

    • w [Numo::NArray] – The eigenvalues.

    • v [Numo::NArray] – The eigenvectors if vals_only is false, otherwise nil.



672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
 
# File 'lib/numo/linalg/function.rb', line 672

def eigh(a, b=nil, vals_only:false, vals_range:nil, uplo:'U', turbo:false)
  jobz = vals_only ? 'N' : 'V' # jobz: Compute eigenvalues and eigenvectors.
  b = a.class.eye(a.shape[0]) if b.nil?
  func = blas_char(a, b) =~ /c|z/ ? 'hegv' : 'sygv'
  if vals_range.nil?
    func << 'd' if turbo
    v, u_, w, = Lapack.call(func.to_sym, a, b, uplo: uplo, jobz: jobz)
    v = nil if vals_only
    [w, v]
  else
    func << 'x'
    il = vals_range.first(1)[0]
    iu = vals_range.last(1)[0]
    a_, b_, w, v, = Lapack.call(func.to_sym, a, b, uplo: uplo, jobz: jobz, range: 'I', il: il + 1, iu: iu + 1)
    v = nil if vals_only
    [w, v]
  end
end

.eigvals(a) ⇒ Numo::NArray

Computes the eigenvalues only for a square nonsymmetric matrix A.

Parameters:

  • a (Numo::NArray) —

    square nonsymmetric matrix (>= 2-dimensinal NArray)

Returns:

  • (Numo::NArray) —

    eigenvalues



696
697
698
699
700
701
702
703
704
705
706
 
# File 'lib/numo/linalg/function.rb', line 696

def eigvals(a)
  jobvl, jobvr = 'N','N'
  case blas_char(a)
  when /c|z/
    w, = Lapack.call(:geev, a, jobvl:jobvl, jobvr:jobvr)
  else
    wr, wi, = Lapack.call(:geev, a, jobvl:jobvl, jobvr:jobvr)
    w  = wr + wi * Complex::I
  end
  w
end

.eigvalsh(a, b = nil, vals_range: nil, uplo: 'U', turbo: false) ⇒ Numo::NArray

Obtains the eigenvalues by solving an ordinary or generalized eigenvalue problem for a square symmetric / Hermitian matrix.

Parameters:

  • a (Numo::NArray) —

    square symmetric/hermitian matrix (>= 2-dimensinal NArray)

  • b (Numo::NArray) (defaults to: nil) —

    (optional) square symmetric matrix (>= 2-dimensinal NArray) If nil, identity matrix is assumed.

  • vals_range (Range) (defaults to: nil) —

    (optional) The range of indices of the eigenvalues (in ascending order) to be returned. If nil or 0…N (N is the size of the matrix a), all eigenvalues are returned.

  • uplo (String or Symbol) (defaults to: 'U') —

    (optional, default=‘U’) Access upper (‘U’) or lower (‘L’) triangle.

Returns:

  • (Numo::NArray) —

    eigenvalues



722
723
724
 
# File 'lib/numo/linalg/function.rb', line 722

def eigvalsh(a, b=nil, vals_range:nil, uplo:'U', turbo:false)
  eigh(a, b, vals_only: true, vals_range: vals_range, uplo: uplo, turbo: turbo).first
end

.expm(a, ord = 8) ⇒ Numo::NArray

Compute the matrix exponential using Pade approximation method.

Examples:

a = Numo::Linalg.expm(Numo::DFloat.zeros([2,2]))
=> Numo::DFloat#shape=[2,2]
[[1, 0],
 [0, 1]]
b = Numo::Linalg.expm(Numo::DFloat[[1, 2], [-1, 3]] * Complex::I)
=> Numo::DComplex#shape=[2,2]
[[0.426459+1.89218i, -2.13721-0.978113i],
 [1.06861+0.489056i, -1.71076+0.914063i]]

Parameters:

  • a (Numo::NArray) —

    square matrix (>= 2-dimensinal NArray)

  • ord (Integer) (defaults to: 8) —

    order of approximation

Returns:

  • (Numo::NArray)

Raises:

  • (NArray::ShapeError) 


1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
 
# File 'lib/numo/linalg/function.rb', line 1192

def expm(a, ord=8)
  raise NArray::ShapeError, 'matrix a is not square matrix' if a.shape[0] != a.shape[1]

  inf_norm = norm(a, 'inf')
  n_squarings = inf_norm.zero? ? 0 : [0, Math.log2(inf_norm).ceil.to_i].max
  a = a / (2**n_squarings)

  sz_mat = a.shape[0]
  c = 1
  s = -1
  x = Numo::DFloat.eye(sz_mat)
  n = Numo::DFloat.eye(sz_mat)
  d = Numo::DFloat.eye(sz_mat)

  (1..ord).each do |k|
    c *= (ord - k + 1).fdiv((2 * ord - k + 1) * k)
    x = a.dot(x)
    cx = c * x
    n += cx
    d += s * cx
    s *= -1
  end

  res = solve(d, n)
  n_squarings.times { res = res.dot(res) }
  res
end

.inv(a, driver: "getrf", uplo: 'U') ⇒ Numo::NArray

Inverse matrix from square matrix ‘a`

Examples:

Numo::Linalg.inv(a,driver:'getrf')
=> Numo::DFloat#shape=[2,2]
[[-2, 1],
 [1.5, -0.5]]
a.dot(Numo::Linalg.inv(a,driver:'getrf'))
=> Numo::DFloat#shape=[2,2]
[[1, 0],
 [8.88178e-16, 1]]

Parameters:

  • a (Numo::NArray) —

    n-by-n square matrix (>= 2-dimensinal NArray)

  • driver (String or Symbol) (defaults to: "getrf") —

    choose LAPACK diriver (‘ge’|‘sy’|‘he’|‘po’) + (“sv”|“trf”) (optional, default=‘getrf’)

  • uplo (String or Symbol) (defaults to: 'U') —

    optional, default=‘U’. Access upper or (‘U’) lower (‘L’) triangle. (omitted when driver:“ge”)

Returns:

  • (Numo::NArray) —

    The inverse matrix.



1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
 
# File 'lib/numo/linalg/function.rb', line 1014

def inv(a, driver:"getrf", uplo:'U')
  case driver
  when /(ge|sy|he|po)sv$/
    d = $1
    b = a.new_zeros.eye
    solve(a, b, driver:d, uplo:uplo)
  when /(ge|sy|he)tr[fi]$/
    d = $1
    lu, piv = lu_fact(a)
    lu_inv(lu, piv)
  when /potr[fi]$/
    lu = cho_fact(a, uplo:uplo)
    cho_inv(lu, uplo:uplo)
  else
    raise ArgumentError, "invalid driver: #{driver}"
  end
end

.ldl(a, uplo: 'U', hermitian: true) ⇒ [lu,d,perm]

Computes the LDLt or Bunch-Kaufman factorization of a symmetric/Hermitian matrix A. The factorization has the form

A = U*D*U**T  or  A = L*D*L**T

where U (or L) is a product of permutation and unit upper (lower) triangular matrices and D is symmetric and block diagonal with 1-by-1 and 2-by-2 diagonal blocks.

Parameters:

  • a (Numo::NArray) —

    m-by-m matrix A (>= 2-dimensinal NArray)

  • uplo (String or Symbol) (defaults to: 'U') —

    optional, default=‘U’. Access upper or (‘U’) lower (‘L’) triangle.

  • hermitian (Bool) (defaults to: true) —

    optional, default=true. If true, hermitian matrix is assumed. (omitted when real-value matrix is given)

Returns:

  • ([lu,d,perm]) —

    • lu [Numo::NArray] – The permutated upper (lower) triangular matrix U (L).

    • d [Numo::NArray] – The block diagonal matrix D.

    • perm [Numo::NArray] – The row-permutation index for changing lu into triangular form.

Raises:

  • (NArray::ShapeError) 


513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
 
# File 'lib/numo/linalg/function.rb', line 513

def ldl(a, uplo: 'U', hermitian: true)
  raise NArray::ShapeError, '2-d array is required' if a.ndim < 2
  raise NArray::ShapeError, 'matrix a is not square matrix' if a.shape[0] != a.shape[1]

  is_complex = blas_char(a) =~ /c|z/
  func = is_complex && hermitian ? 'hetrf' : 'sytrf'
  lud, ipiv, = Lapack.call(func.to_sym, a, uplo: uplo)

  lu = (uplo == 'U' ? lud.triu : lud.tril).tap { |mat| mat[mat.diag_indices(0)] = 1.0 }
  d = lud[lud.diag_indices(0)].diag

  m = a.shape[0]
  n = m - 1
  changed_2x2 = false
  perm = Numo::Int32.new(m).seq
  m.times do |t|
    i = uplo == 'U' ? t : n - t
    j = uplo == 'U' ? i - 1 : i + 1;
    r = uplo == 'U' ? 0..i : i..n;
    if ipiv[i] > 0
      k = ipiv[i] - 1
      lu[[k, i], r] = lu[[i, k], r].dup
      perm[[k, i]] = perm[[i, k]].dup
    elsif j.between?(0, n) && ipiv[i] == ipiv[j] && !changed_2x2
      k = ipiv[i].abs - 1
      d[j, i] = lud[j, i]
      d[i, j] = is_complex && hermitian ? lud[j, i].conj : lud[j, i]
      lu[j, i] = 0.0
      lu[[k, j], r] = lu[[j, k], r].dup
      perm[[k, j]] = perm[[j, k]].dup
      changed_2x2 = true
      next
    end
    changed_2x2 = false
  end

  [lu, d, perm.sort_index]
end

.lstsq(a, b, driver: 'lsd', rcond: -1)) ⇒ [x, resids, rank, s]

Computes the minimum-norm solution to a linear least squares problem:

minimize 2-norm(| b - A*x |)

using the singular value decomposition (SVD) of A. A is an M-by-N matrix which may be rank-deficient.

Parameters:

  • a (Numo::NArray) —

    m-by-n matrix A (>= 2-dimensinal NArray)

  • b (Numo::NArray) —

    m-by-nrhs right-hand-side matrix b (>= 1-dimensinal NArray)

  • driver (String or Symbol) (defaults to: 'lsd') —

    choose LAPACK driver from ‘lsd’,‘lss’,‘lsy’ (optional, default=‘lsd’)

  • rcond (Float) (defaults to: -1)) —

    (optional, default=-1) RCOND is used to determine the effective rank of A. Singular values ‘S(i) <= RCOND*S(1)` are treated as zero. If RCOND < 0, machine precision is used instead.

Returns:

  • ([x, resids, rank, s]) —

    • x – The solution matrix/vector X.

    • resids – Sums of residues, squared 2-norm for each column in ‘b - a x`. If matrix_rank(a) < N or > M, or ’gelsy’ is used, this is an empty array.

    • rank – The effective rank of A, i.e., the number of singular values which are greater than RCOND*S(1).

    • s – The singular values of A in decreasing order. Returns nil if ‘gelsy’ is used.



1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
 
# File 'lib/numo/linalg/function.rb', line 1058

def lstsq(a, b, driver:'lsd', rcond:-1)
  a = NArray.asarray(a)
  b = NArray.asarray(b)
  b_orig = nil
  if b.shape.size==1
    b_orig = b
    b = b_orig[true,:new]
  end
  m = a.shape[-2]
  n = a.shape[-1]
  #nrhs = b.shape[-1]
  if m != b.shape[-2]
    raise NArray::ShapeError, "size mismatch: A-row and B-row"
  end
  if m < n   # need to extend b matrix
    shp = b.shape
    shp[-2] = n
    b2 = b.class.zeros(*shp)
    b2[false,0...m,true] = b
    b = b2
  end
  case driver.to_s
  when /^(ge)?lsd$/i
    # x, s, rank, info
    x, s, rank, = Lapack.call(:gelsd, a, b, rcond:rcond)
  when /^(ge)?lss$/i
    # v, x, s, rank, info
    _, x, s, rank, = Lapack.call(:gelss, a, b, rcond:rcond)
  when /^(ge)?lsy$/i
    jpvt = Int32.zeros(*a[false,0,true].shape)
    # v, x, jpvt, rank, info
    _, x, _, rank, = Lapack.call(:gelsy, a, b, jpvt, rcond:rcond)
    s = nil
  else
    raise ArgumentError, "invalid driver: #{driver}"
  end
  resids = nil
  if m > n
    if /ls(d|s)$/i =~ driver
      case rank
      when n
        resids = (x[n..-1,true].abs**2).sum(axis:0)
      when NArray
        resids = (x[false,n..-1,true].abs**2).sum(axis:-2)
        ## NArray does not suppurt this yet.
        # resids = x[false,0,true].new_zeros
        # mask = rank.eq(n)
        # resids[mask,true] = (x[mask,n..-1,true].abs**2).sum(axis:-2)
      end
    end
    x = x[false,0...n,true]
  end
  if b_orig && b_orig.shape.size==1
    x = x[true,0]
    resids &&= resids[false,0]
  end
  [x, resids, rank, s]
end

.lu(a, permute_l: false) ⇒ [p,l,u], [pl,u]

Computes an LU factorization of a M-by-N matrix A using partial pivoting with row interchanges.

The factorization has the form

A = P * L * U

where P is a permutation matrix, L is lower triangular with unit diagonal elements (lower trapezoidal if m > n), and U is upper triangular (upper trapezoidal if m < n).

Parameters:

  • a (Numo::NArray) —

    m-by-n matrix A (>= 2-dimensinal NArray)

  • permute_l (Bool) (defaults to: false) —

    (optional) If true, perform the matrix product of P and L.

Returns:

  • ([p,l,u]) —

    if permute_l == false

  • ([pl,u]) —

    if permute_l == true

    • p [Numo::NArray] – The permutation matrix P.

    • l [Numo::NArray] – The factor L.

    • u [Numo::NArray] – The factor U.

Raises:

  • (NArray::ShapeError) 


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

def lu(a, permute_l: false)
  raise NArray::ShapeError, '2-d array is required' if a.ndim < 2
  m, n = a.shape
  k = [m, n].min
  lu, ip = lu_fact(a)
  l = lu.tril.tap { |mat| mat[mat.diag_indices(0)] = 1.0 }[true, 0...k]
  u = lu.triu[0...k, 0...n]
  p = Numo::DFloat.eye(m).tap do |mat|
        ip.to_a.each_with_index { |i, j| mat[true, [i - 1, j]] = mat[true, [j, i - 1]].dup }
      end
  permute_l ? [p.dot(l), u] : [p, l, u]
end

.lu_fact(a) ⇒ [lu, ipiv]

Computes an LU factorization of a M-by-N matrix A using partial pivoting with row interchanges.

The factorization has the form

A = P * L * U

where P is a permutation matrix, L is lower triangular with unit diagonal elements (lower trapezoidal if m > n), and U is upper triangular (upper trapezoidal if m < n).

Parameters:

  • a (Numo::NArray) —

    m-by-n matrix A (>= 2-dimensinal NArray)

Returns:

  • ([lu, ipiv]) —

    • lu [Numo::NArray] – The factors L and U from the factorization ‘A = P*L*U`; the unit diagonal elements of L are not stored.

    • ipiv [Numo::NArray] – The pivot indices; for 1 <= i <= min(M,N),

      row i of the matrix was interchanged with row IPIV(i).



444
445
446
 
# File 'lib/numo/linalg/function.rb', line 444

def lu_fact(a)
  Lapack.call(:getrf, a)[0..1]
end

.lu_inv(lu, ipiv) ⇒ Numo::NArray

Computes the inverse of a matrix using the LU factorization computed by Numo::Linalg.lu_fact.

This method inverts U and then computes inv(A) by solving the system

inv(A)*L = inv(U)

for inv(A).

Parameters:

  • lu (Numo::NArray) —

    matrix containing the factors L and U from the factorization ‘A = P*L*U` as computed by Numo::Linalg.lu_fact.

  • ipiv (Numo::NArray) —

    The pivot indices from Numo::Linalg.lu_fact; for 1<=i<=N, row i of the matrix was interchanged with row IPIV(i).

Returns:

  • (Numo::NArray) —

    the inverse of the original matrix A.



465
466
467
 
# File 'lib/numo/linalg/function.rb', line 465

def lu_inv(lu, ipiv)
  Lapack.call(:getri, lu, ipiv)[0]
end

.lu_solve(lu, ipiv, b, trans: "N") ⇒ Numo::NArray

Solves a system of linear equations

A * X = B  or  A**T * X = B

with a N-by-N matrix A using the LU factorization computed by Numo::Linalg.lu_fact

Parameters:

  • lu (Numo::NArray) —

    matrix containing the factors L and U from the factorization ‘A = P*L*U` as computed by Numo::Linalg.lu_fact.

  • ipiv (Numo::NArray) —

    The pivot indices from Numo::Linalg.lu_fact; for 1<=i<=N, row i of the matrix was interchanged with row IPIV(i).

  • b (Numo::NArray) —

    the right hand side matrix B.

  • trans (String or Symbol) (defaults to: "N") —

    Specifies the form of the system of equations:

    - If 'N': `A * X = B` (No transpose).
- If 'T': `A*\*T* X = B` (Transpose).
- If 'C': `A*\*T* X = B` (Conjugate transpose = Transpose).

Returns:

  • (Numo::NArray) —

    the solution matrix X.



490
491
492
 
# File 'lib/numo/linalg/function.rb', line 490

def lu_solve(lu, ipiv, b, trans:"N")
  Lapack.call(:getrs, lu, ipiv, b, trans:trans)[0]
end

.matmul(a, b) ⇒ Numo::NArray

Matrix product.

Parameters:

  • a (Numo::NArray) —

    matrix (>= 2-dimensinal NArray)

  • b (Numo::NArray) —

    matrix (>= 2-dimensinal NArray)

Returns:

  • (Numo::NArray) —

    result of matrix product



166
167
168
 
# File 'lib/numo/linalg/function.rb', line 166

def matmul(a, b)
  Blas.call(:gemm, a, b)
end

.matrix_power(a, n) ⇒ Object

Compute a square matrix ‘a` to the power `n`.

* If n > 0: return `a**n`.
* If n == 0: return identity matrix.
* If n < 0: return `(a*\*-1)*\*n.abs`.

Examples:

i = Numo::DFloat[[0, 1], [-1, 0]]
=> Numo::DFloat#shape=[2,2]
[[0, 1],
 [-1, 0]]
Numo::Linalg.matrix_power(i,3)
=> Numo::DFloat#shape=[2,2]
[[0, -1],
 [1, 0]]
Numo::Linalg.matrix_power(i,0)
=> Numo::DFloat#shape=[2,2]
[[1, 0],
 [0, 1]]
Numo::Linalg.matrix_power(i,-3)
=> Numo::DFloat#shape=[2,2]
[[0, 1],
 [-1, 0]]

q = Numo::DFloat.zeros(4,4)
q[0..1,0..1] = -i
q[2..3,2..3] = i
q
=> Numo::DFloat#shape=[4,4]
[[-0, -1, 0, 0],
 [1, -0, 0, 0],
 [0, 0, 0, 1],
 [0, 0, -1, 0]]
Numo::Linalg.matrix_power(q,2)
=> Numo::DFloat#shape=[4,4]
[[-1, 0, 0, 0],
 [0, -1, 0, 0],
 [0, 0, -1, 0],
 [0, 0, 0, -1]]

Parameters:

  • a (Numo::NArray) —

    square matrix (>= 2-dimensinal NArray).

  • n (Integer) —

    the exponent.



212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
 
# File 'lib/numo/linalg/function.rb', line 212

def matrix_power(a, n)
  a = NArray.asarray(a)
  m,k = a.shape[-2..-1]
  unless m==k
    raise NArray::ShapeError, "input must be a square array"
  end
  unless Integer===n
    raise ArgumentError, "exponent must be an integer"
  end
  if n == 0
    return a.class.eye(m)
  elsif n < 0
    a = inv(a)
    n = n.abs
  end
  if n <= 3
    r = a
    (n-1).times do
      r = matmul(r,a)
    end
  else
    while (n & 1) == 0
      a = matmul(a,a)
      n >>= 1
    end
    r = a
    while n != 0
      a = matmul(a,a)
      n >>= 1
      if (n & 1) != 0
        r = matmul(r,a)
      end
    end
  end
  r
end

.matrix_rank(m, tol: nil, driver: 'svd') ⇒ Object

Compute matrix rank of array using SVD Rank is the number of singular values greater than tol.

Parameters:

  • m (Numo::NArray) —

    matrix (>= 2-dimensional NArray)

  • tol (Float) (defaults to: nil) —

    threshold below which singular values are considered to be zero. If tol is nil, ‘tol = sing_vals.max() * m.shape.max * EPSILON`.

  • driver (String or Symbol) (defaults to: 'svd') —

    choose LAPACK solver from ‘svd’, ‘sdd’. (optional, default=‘svd’)



950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
 
# File 'lib/numo/linalg/function.rb', line 950

def matrix_rank(m, tol:nil, driver:'svd')
  m = Numo::NArray.asarray(m)
  if m.ndim < 2
    m.ne(0).any? ? 1 : 0
  else
    case driver.to_s
    when /^(ge)?sdd$/, "turbo"
      s = Lapack.call(:gesdd, m, jobz:'N')[0]
    when /^(ge)?svd$/
      s = Lapack.call(:gesvd, m, jobu:'N', jobvt:'N')[0]
    else
      raise ArgumentError, "invalid driver: #{driver}"
    end
    tol ||= s.max(axis:-1, keepdims:true) *
      (m.shape[-2..-1].max * s.class::EPSILON)
    (s > tol).count(axis:-1)
  end
end

.norm(a, ord = nil, axis: nil, keepdims: false) ⇒ Numo::NArray

Compute matrix or vector norm.

|  ord  |  matrix norm           | vector norm                 |
| ----- | ---------------------- | --------------------------- |
|  nil  | Frobenius norm         | 2-norm                      |
| 'fro' | Frobenius norm         |  -                          |
| 'inf' | x.abs.sum(axis:-1).max | x.abs.max                   |
|    0  |  -                     | (x.ne 0).sum                |
|    1  | x.abs.sum(axis:-2).max | same as below               |
|    2  | 2-norm (max sing_vals) | same as below               |
| other |  -                     | (x.abs**ord).sum**(1.0/ord) |

Parameters:

  • a (Numo::NArray) —

    matrix or vector (>= 1-dimensinal NArray)

  • ord (String or Symbol) (defaults to: nil) —

    Order of the norm .

  • axis (Integer or Array) (defaults to: nil) —

    Applied axes (optional).

  • keepdims (Bool) (defaults to: false) —

    If true, the applied axes are left in result with size one (optional).

Returns:

  • (Numo::NArray) —

    norm result



748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
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
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
 
# File 'lib/numo/linalg/function.rb', line 748

def norm(a, ord=nil, axis:nil, keepdims:false)
  a = Numo::NArray.asarray(a)

  # check axis
  if axis
    case axis
    when Integer
      axis = [axis]
    when Array
      if axis.size < 1 || axis.size > 2
        raise ArgmentError, "axis option should be 1- or 2-element array"
      end
    else
      raise ArgumentError, "invalid option for axis: #{axis}"
    end
    # swap axes
    if a.ndim > 1
      idx = (0...a.ndim).to_a
      tmp = []
      axis.each do |i|
        x = idx[i]
        if x.nil?
          raise ArgmentError, "axis contains same dimension"
        end
        tmp << x
        idx[i] = nil
      end
      idx.compact!
      idx.concat(tmp)
      a = a.transpose(*idx)
    end
  else
    case a.ndim
    when 0
      raise ArgumentError, "zero-dimensional array"
    when 1
      axis = [-1]
    else
      axis = [-2,-1]
    end
  end

  # calculate norm
  case axis.size

  when 1  # vector
    k = keepdims
    ord ||= 2  # default
    case ord.to_s
    when "0"
      r = a.class.cast(a.ne(0)).sum(axis:-1, keepdims:k)
    when "1"
      r = a.abs.sum(axis:-1, keepdims:k)
    when "2"
      r = Blas.call(:nrm2, a, keepdims:k)
    when /^-?\d+$/
      o = ord.to_i
      r = (a.abs**o).sum(axis:-1, keepdims:k)**(1.0/o)
    when /^inf(inity)?$/i
      r = a.abs.max(axis:-1, keepdims:k)
    when /^-inf(inity)?$/i
      r = a.abs.min(axis:-1, keepdims:k)
    else
      raise ArgumentError, "ord (#{ord}) is invalid for vector norm"
    end

  when 2  # matrix
    if keepdims
      fixdims = [true] * a.ndim
      axis.each do |i|
        if i < -a.ndim || i >= a.ndim
          raise ArgmentError, "axis (%d) is out of range", i
        end
        fixdims[i] = :new
      end
    end
    ord ||= "fro"  # default
    case ord.to_s
    when "1"
      r, = Lapack.call(:lange, a, '1')
    when "-1"
      r = a.abs.sum(axis:-2).min(axis:-1)
    when "2"
      svd, = Lapack.call(:gesvd, a, jobu:'N', jobvt:'N')
      r = svd.max(axis:-1)
    when "-2"
      svd, = Lapack.call(:gesvd, a, jobu:'N', jobvt:'N')
      r = svd.min(axis:-1)
    when /^f(ro)?$/i
      r, = Lapack.call(:lange, a, 'F')
    when /^inf(inity)?$/i
      r, = Lapack.call(:lange, a, 'I')
    when /^-inf(inity)?$/i
      r = a.abs.sum(axis:-1).min(axis:-1)
    else
      raise ArgumentError, "ord (#{ord}) is invalid for matrix norm"
    end
    if keepdims
      if NArray===r
        r = r[*fixdims]
      else
        r = a.class.new(1,1).store(r)
      end
    end
  end
  return r
end

.null_space(a, rcond: -1)) ⇒ Numo::NArray

Computes an orthonormal basis for the null space of matrix A.

Parameters:

  • a (Numo::NArray) —

    m-by-n matrix A (>= 2-dimensional NArray).

  • rcond (Float) (defaults to: -1)) —

    (optional) rcond is used to determine the effective rank of A. Singular values ‘s <= rcond * s.max` are treated as zero. If rcond < 0, machine precision is used instead.

Returns:

  • (Numo::NArray) —

    The orthonormal basis for the null space of matrix A.

Raises:

  • (NArray::ShapeError) 


383
384
385
386
387
388
389
390
391
 
# File 'lib/numo/linalg/function.rb', line 383

def null_space(a, rcond: -1)
  raise NArray::ShapeError, '2-d array is required' if a.ndim < 2
  s, _u, vh = svd(a)
  tol = s.max * (rcond.nil? || rcond < 0 ? a.class::EPSILON * a.shape.max : rcond)
  k = (s > tol).count
  return a.class.new if k == vh.shape[0]
  r = vh[k..-1, true].transpose.dup
  blas_char(vh) =~ /c|z/ ? r.conj : r
end

.orth(a, rcond: -1)) ⇒ Numo::NArray

Computes an orthonormal basis for the range of matrix A.

Parameters:

  • a (Numo::NArray) —

    m-by-n matrix A (>= 2-dimensional NArray).

  • rcond (Float) (defaults to: -1)) —

    (optional) rcond is used to determine the effective rank of A. Singular values ‘s <= rcond * s.max` are treated as zero. If rcond < 0, machine precision is used instead.

Returns:

  • (Numo::NArray) —

    The orthonormal basis for the range of matrix A.

Raises:

  • (NArray::ShapeError) 


366
367
368
369
370
371
372
 
# File 'lib/numo/linalg/function.rb', line 366

def orth(a, rcond: -1)
  raise NArray::ShapeError, '2-d array is required' if a.ndim < 2
  s, u, = svd(a)
  tol = s.max * (rcond.nil? || rcond < 0 ? a.class::EPSILON * a.shape.max : rcond)
  k = (s > tol).count
  u[true, 0...k]
end

.pinv(a, driver: "svd", rcond: nil) ⇒ Numo::NArray

Compute the (Moore-Penrose) pseudo-inverse of a matrix using svd or lstsq.

Examples:

a = Numo::DFloat.new(5,3).rand_norm
=> Numo::DFloat#shape=[5,3]
[[-0.581255, -0.168354, 0.586895],
 [-0.595142, -0.802802, -0.326106],
 [0.282922, 1.68427, 0.918499],
 [-0.0485384, -0.464453, -0.992194],
 [0.413794, -0.60717, -0.699695]]
b = Numo::Linalg.pinv(a,driver:"svd")
=> Numo::DFloat(view)#shape=[3,5]
[[-0.360863, -0.813125, -0.353367, -0.891963, 0.877253],
 [-0.227645, 0.162939, 0.696655, 0.787685, -0.469346],
 [0.408671, -0.308323, -0.337807, -1.13833, 0.228051]]
(a-a.dot(b.dot(a))).abs.max
=> 5.551115123125783e-16

Parameters:

  • a (Numo::NArray) —

    m-by-n matrix A (>= 2-dimensinal NArray)

  • driver (String or Symbol) (defaults to: "svd") —

    choose LAPACK driver from SVD (‘svd’, ‘sdd’) or Least square (‘lsd’,‘lss’,‘lsy’) (optional, default=‘svd’)

  • rcond (Float) (defaults to: nil) —

    (optional, default=-1) RCOND is used to determine the effective rank of A. Singular values ‘S(i) <= RCOND*S(1)` are treated as zero. If RCOND < 0, machine precision is used instead.

Returns:

  • (Numo::NArray) 


1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
 
# File 'lib/numo/linalg/function.rb', line 1145

def pinv(a, driver:"svd", rcond:nil)
  a = NArray.asarray(a)
  if a.ndim < 2
    raise NArray::ShapeError, "2-d array is required"
  end
  case driver
  when /^(ge)?s[dv]d$/
    s, u, vh = svd(a, driver:driver, job:'S')
    if rcond.nil? || rcond < 0
      rcond = ((SFloat===s) ? 1e3 : 1e6) * s.class::EPSILON
    elsif ! Numeric === rcond
      raise ArgumentError, "rcond must be Numeric"
    end
    cond = (s > rcond * s.max(axis:-1, keepdims:true))
    if cond.all?
      r = s.reciprocal
    else
      r = s.new_zeros
      r[cond] = s[cond].reciprocal
    end
    u *= r[false,:new,true]
    dot(u,vh).conj.swapaxes(-2,-1)
  when /^(ge)?ls[dsy]$/
    b = a.class.eye(a.shape[-2])
    x, = lstsq(a, b, driver:driver, rcond:rcond)
    x
  else
    raise ArgumentError, "#{driver.inspect} is not one of drivers: "+
      "svd, sdd, lsd, lss, lsy"
  end
end

.qr(a, mode: "reduce") ⇒ r, ...

Computes a QR factorization of a complex M-by-N matrix A: A = Q * R.

Parameters:

  • a (Numo::NArray) —

    m-by-n matrix A (>= 2-dimensinal NArray)

  • mode (String) (defaults to: "reduce") —

    • “reduce” – returns both Q and R,

    • “r” – returns only R,

    • “economic” – returns both Q and R but computed in economy-size,

    • “raw” – returns QR and TAU used in LAPACK.

Returns:

  • (r) —

    if mode:“r”

  • ([q,r]) —

    if mode:“reduce” or “economic”

  • ([qr,tau]) —

    if mode:“raw” (LAPACK geqrf result)



264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
 
# File 'lib/numo/linalg/function.rb', line 264

def qr(a, mode:"reduce")
  qr,tau, = Lapack.call(:geqrf, a)
  *shp,m,n = qr.shape
  r = (m >= n && %w[economic raw].include?(mode)) ?
    qr[false, 0...n, true].triu : qr.triu
  mode = mode.to_s.downcase
  case mode
  when "r"
    return r
  when "raw"
    return [qr,tau]
  when "reduce","economic"
    # skip
  else
    raise ArgumentError, "invalid mode:#{mode}"
  end
  if m < n
    q, = Lapack.call(:orgqr, qr[false, 0...m], tau)
  elsif mode == "economic"
    q, = Lapack.call(:orgqr, qr, tau)
  else
    qqr = qr.class.zeros(*(shp+[m,m]))
    qqr[false,0...n] = qr
    q, = Lapack.call(:orgqr, qqr, tau)
  end
  return [q,r]
end

.slogdet(a) ⇒ [sign,logdet]

Natural logarithm of the determinant of a matrix

Parameters:

  • a (Numo::NArray) —

    matrix (>= 2-dimensional NArray)

Returns:

  • ([sign,logdet]) —

    • sign – A number representing the sign of the determinant.

    • logdet – The natural log of the absolute value of the determinant.



925
926
927
928
929
930
931
932
933
934
935
936
937
938
 
# File 'lib/numo/linalg/function.rb', line 925

def slogdet(a)
  lu, piv, = Lapack.call(:getrf, a)
  idx = piv.new_narray.store(piv.class.new(piv.shape[-1]).seq(1))
  m = piv.eq(idx).count_false(axis:-1) % 2
  sign = m * -2 + 1

  lud = lu.diagonal
  if (lud.eq 0).any?
    return 0, (-Float::INFINITY)
  end
  lud_abs = lud.abs
  sign *= (lud/lud_abs).prod
  [sign, NMath.log(lud_abs).sum(axis:-1)]
end

.solve(a, b, driver: "gen", uplo: 'U') ⇒ Numo::NArray

Solves linear equation ‘a * x = b` for `x` from square matrix `a`

Parameters:

  • a (Numo::NArray) —

    n-by-n square matrix (>= 2-dimensinal NArray)

  • b (Numo::NArray) —

    n-by-nrhs right-hand-side matrix (>= 1-dimensinal NArray)

  • driver (String or Symbol) (defaults to: "gen") —

    choose LAPACK diriver from ‘gen’,‘sym’,‘her’ or ‘pos’. (optional, default=‘gen’)

  • uplo (String or Symbol) (defaults to: 'U') —

    optional, default=‘U’. Access upper or (‘U’) lower (‘L’) triangle. (omitted when driver:“gen”)

Returns:

  • (Numo::NArray) —

    The solusion matrix/vector X.



983
984
985
986
987
988
989
990
991
992
993
994
 
# File 'lib/numo/linalg/function.rb', line 983

def solve(a, b, driver:"gen", uplo:'U')
  case driver.to_s
  when /^gen?(sv)?$/i
    # returns lu, x, ipiv, info
    Lapack.call(:gesv, a, b)[1]
  when /^(sym?|her?|pos?)(sv)?$/i
    func = driver[0..1].downcase+"sv"
    Lapack.call(func, a, b, uplo:uplo)[1]
  else
    raise ArgumentError, "invalid driver: #{driver}"
  end
end

.svd(a, driver: 'svd', job: 'A') ⇒ [sigma,u,vt]

Computes the Singular Value Decomposition (SVD) of a M-by-N matrix A, and the left and/or right singular vectors. The SVD is written

A = U * SIGMA * transpose(V)

where SIGMA is an M-by-N matrix which is zero except for its min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA are the singular values of A; they are real and non-negative, and are returned in descending order. The first min(m,n) columns of U and V are the left and right singular vectors of A. Note that the routine returns V**T, not V.

Parameters:

  • a (Numo::NArray) —

    m-by-n matrix A (>= 2-dimensinal NArray)

  • driver (String or Symbol) (defaults to: 'svd') —

    choose LAPACK solver from ‘svd’, ‘sdd’. (optional, default=‘svd’)

  • job (String or Symbol) (defaults to: 'A') —

    • ‘A’: all M columns of U and all N rows of V**T are returned in the arrays U and VT.

    • ‘S’: the first min(M,N) columns of U and the first min(M,N) rows of V**T are returned in the arrays U and VT.

    • ‘N’: no columns of U or rows of V**T are computed.

Returns:

  • ([sigma,u,vt]) —

    SVD result. Array<Numo::NArray>



317
318
319
320
321
322
323
324
325
326
327
328
329
 
# File 'lib/numo/linalg/function.rb', line 317

def svd(a, driver:'svd', job:'A')
  unless /^[ASN]/i =~ job
    raise ArgumentError, "invalid job: #{job.inspect}"
  end
  case driver.to_s
  when /^(ge)?sdd$/i, "turbo"
    Lapack.call(:gesdd, a, jobz:job)[0..2]
  when /^(ge)?svd$/i
    Lapack.call(:gesvd, a, jobu:job, jobvt:job)[0..2]
  else
    raise ArgumentError, "invalid driver: #{driver}"
  end
end

.svdvals(a, driver: 'svd') ⇒ Numo::NArray

Computes the Singular Values of a M-by-N matrix A. The SVD is written

A = U * SIGMA * transpose(V)

where SIGMA is an M-by-N matrix which is zero except for its min(m,n) diagonal elements. The diagonal elements of SIGMA are the singular values of A; they are real and non-negative, and are returned in descending order.

Parameters:

  • a (Numo::NArray) —

    m-by-n matrix A (>= 2-dimensinal NArray)

  • driver (String or Symbol) (defaults to: 'svd') —

    choose LAPACK solver from ‘svd’, ‘sdd’. (optional, default=‘svd’)

Returns:

  • (Numo::NArray) —

    returns SIGMA (singular values).



346
347
348
349
350
351
352
353
354
355
 
# File 'lib/numo/linalg/function.rb', line 346

def svdvals(a, driver:'svd')
  case driver.to_s
  when /^(ge)?sdd$/i, "turbo"
    Lapack.call(:gesdd, a, jobz:'N')[0]
  when /^(ge)?svd$/i
    Lapack.call(:gesvd, a, jobu:'N', jobvt:'N')[0]
  else
    raise ArgumentError, "invalid driver: #{driver}"
  end
end