Linear Algebra

Module/Classes

GSL (Module)
- Linalg (Module)
  - LU (Module)
    - LUMatrix < Matrix
  - QR (Module)
    - QRMatrix < Matrix
    - QMatrix < Matrix
    - RMatrix < Matrix
    - TauVector < Vector
  - QRPT (Module)
    - QRMatrix < Matrix
    - QMatrix < Matrix
    - RMatrix < Matrix
    - TauVector < Vector
  - SV (Module)
    - UMatrix < Matrix
    - VMatrix < Matrix
    - SVector < Vector
  - Cholesky (Module)
    - CholeskyMatrix < Matrix
  - HH (Module)

LU Decomposition

GSL::Linalg::LU.decomp(A)

GSL::Matrix#LU_decomp

These method calculate the LU decomposition of the matrix. The returned value is an array of [LU, perm, sign].

ex1:

include GSL::Linalg
m = Matrix.new(...)
lu, perm, sign = LU.decomp(m)

ex2:
```
lu, perm, sign = m.LU_decomp
```

GSL::Linalg::LU.solve(A, b)

GSL::Linalg::LU.solve(lu, perm, b)

GSL::Matrix#LU_solve(b)

GSL::Linalg::LU::LUMatrix#solve(perm, b)

The following is an example to solve a linear system

A x = b, b = [1, 2, 3, 4]

using LU decomposition.

ex1)

require 'gsl'
A = Matrix.new([0.18, 0.60, 0.57, 0.96], [0.41, 0.24, 0.99, 0.58],
                   [0.14, 0.30, 0.97, 0.66], [0.51, 0.13, 0.19, 0.85])
lu, perm, sign = A.LU_decomp         
p lu.class                         <- GSL::Linalg::LU::LUMatrix, subclass of the Matrix class
b = [1, 2, 3, 4].to_gv             <- "b" must be given as a GSL::Vector
x = Linalg::LU.solve(lu, perm, b)    
p x                                <- The solution is returned by a vector.

ex2)

A = Matrix.new([0.18, 0.60, 0.57, 0.96], [0.41, 0.24, 0.99, 0.58],
                   [0.14, 0.30, 0.97, 0.66], [0.51, 0.13, 0.19, 0.85])
p A.class                           <- GSL::Matrix
x = Linalg::LU.solve(A, b)          <- LU decomposition is calculated internally (A is not modified)

GSL::Linalg::LU.svx(A, b)

GSL::Linalg::LU.svx(lu, perm, b)

GSL::Matrix#svx(b)

GSL::Linalg::LU::LUMatrix#svx(perm, b)

These solve the system Ax = b. The input vector b is overwitten by the solution x.

GSL::Linalg::LU.refine(A, lu, perm, b, x)

This method applys an iterative improvement to x, the solution of A x = b, using the LU decomposition of A.

GSL::Linalg::LU.invert(A)

GSL::Linalg::LU.invert(lu, perm)

GSL::Matrix#invert

GSL::Linalg::LU::LUMatrix#invert(perm)

These computes and returns the inverse of the matrix.

GSL::Linalg::LU.det(A)

GSL::Linalg::LU.det(lu, signum)

GSL::Matrix#det

GSL::Linalg::LU::LUMatrix#det(signum)

These methods return the determinant of the matrix.

QR decomposition

GSL::Linalg::QR.decomp(A)

GSL::Matrix#QR_decomp

These compute QR decomposition of the matrix and return an array [QR, tau].

ex:

include GSL::Linalg
m = Matrix.new(...)
qr, tau = m.QR_decomp
p qr.class                 <----- GSL::Linalg::QR::QRMatrix, subclass of GSL::Matrix
p tau.class                <----- GSK::Linalg::QR::TauVector, subclass of GSL::Vector

qr, tau = QR.decomp(m)

GSL::Linalg::QR.solve(A, b)

GSL::Linalg::QR.solve(QR, tau, b)

GSL::Matrix#QR_solve(b)

GSL::Linalg::QR::QRMatrix#solve(tau, b)

Solve the system A x = b using the QR decomposition.

ex1:

include GSL::Linalg
m = Matrix.new(...)
b = Vector.new(...)
x = QR.solve(m, b)

ex2:

qr, tau = QR.decomp(m)
x = QR.solve(qr, tau, b)

ex3:
```
x = m.QR_solve(b)
```
ex4:
```
x = qr.solve(tau, b)
```

GSL::Linalg::QR.svx(A, x)

GSL::Linalg::QR.svx(QR, tau, x)

GSL::Matrix#QR_svx(x)

GSL::Linalg::QR::QRMatrix#svx(tau, x)

Solve the system A x = b. The input vector x is first give by the right-hand side vector b, and is overwritten by the solution.

GSL::Linalg::QR.unpack(QR, tau)

GSL::Linalg::QR::QRMatrix#unpack(tau)

Unpack the encoded QR decomposition QR,tau and return an array [Q, R].

GSL::Linalg::QR.QRsolve(Q, R, tau)

This method solves the system R x = Q^T b for x. It can be used when the QR decomposition of a matrix is available in unpacked form as Q,R.

QR Decomposition with Column Pivoting

GSL::Linalg::QRPT.decomp(A)

GSL::Matrix#QRPT_decomp

These methods factorize the M-by-N matrix A into the QRP^T decomposition A = Q R P^T, and return an array [QR, tau, perm, signum].

ex1:

require("gsl")
include GSL::Linalg
m = Matrix.new(...)
qr, tau, perm = QRPT.decomp(m)
p qr.class                 <----- GSL::Linalg::QRPT::QRMatrix, subclass of GSL::Matrix

ex2:
```
qr, tau, perm = m.QROT_decomp
```

GSL::Linalg::QRPT.decomp2(A)

GSL::Matrix#QRPT_decomp2

These return an array [Q, R, tau, perm, signum].

q, r, tau, perm = QRPT.decomp2(m)
p q.class                  <----- GSL::Linalg::QRPT::QMatrix
p r.class                  <----- GSL::Linalg::QRPT::RMatrix

GSL::Linalg::QRPT.solve(m, b)

GSL::Linalg::QRPT.solve(qr, tau, perm, b)

GSL::Matrix#QRPT_solve(A, b)

GSL::Linalg::QRPT::QRMatrix#solve(qr, tau, perm, b)

These methods solve the system A x = b using the QRP^T decomposition of A into QR, tau, perm. The solution x is returned as a Vector.

ex1:

include GSL::Linalg
m = Matrix.new(...)
qr, tau, perm = m.QRPT_decomp
b = Vector.new([1, 2, 3, 4])
x = QRPT.solve(qr, tau, perm, b)

ex2:
```
x = QRPT.solve(m, b)
```
ex3:
```
x = qr.solve(tau, p, b)
```
ex4:
```
x = m.QRPT_solve(b)
```

GSL::Linalg::QRPT.svx(m, b)

GSL::Linalg::QRPT.svx(qr, tau, perm, b)

GSL::Matrix#QRPT_svx(A, b)

These methods solve the system A x = b using the QRP^T decomposition of A into QR, tau, perm. The input b is overwritten by the solution x.

GSL::Linalg::QRPT.QRsolve(q, r, tau, perm, b)

This method solves the system R P^T x = Q^T b for x. It can be used when the QR decomposition of a matrix is available in unpacked form as q, r obtained by the method decomp2.

ex:

q, r, tau, perm = QRPT_decomp2
x = QRPT.QRsolve(q, r, perm, b)

GSL::Linalg::QRPT.update(q, r, perm, u, v)

GSL::Linalg::QRPT.Rsolve(qr, perm, b)

GSL::Linalg::QRPT::QRMatrix#Rsolve(perm, b)

GSL::Linalg::QRPT.Rsvx(qr, perm, b)

GSL::Linalg::QRPT::QRMatrix#Rsvx(perm, b)

Singular Value Decomposition

GSL::Linalg::SV.decomp(A)

GSL::Matrix#SV_decomp

These methods factorize the M-by-N matrix A into the singular value decomposition A = U S V^T using the Golub-Reinsch SVD algorithm, and return an array [U, V, S].

ex1:

include GSL::Linalg
m = Matrix.new(...)
u, v, s = SV.decomp(m)
p u.class              <----- GSL::Linalg::SV::UMatrix
p v.class              <----- GSL::Linalg::SV::VMatrix
p s.class              <----- GSL::Linalg::SV::SVector

ex2:
```
u, v, s = m.SV_decomp
```

GSL::Linalg::SV.decomp_mod(A)

GSL::Matrix#SV_decomp_mod

These compute the SVD using the modified Golub-Reinsch algorithm, which is faster for M>>N.

GSL::Linalg::SV.decomp_jacobi(A)

GSL::Matrix#SV_decomp_jacobi

These compute the SVD using one-sided Jacobi orthogonalization. The Jacobi method can compute singular values to higher relative accuracy than Golub-Reinsch algorithms.

GSL::Linalg::SV.solve(A, b)

GSL::Linalg::SV.solve(U, V, S, b)

GSL::Matrix#SV_solve(b)

These methods solve the system A x = b using the singular value decomposition U, S, V of A.

ex1:

include GSL::Linalg
m = Matrix.new(...)
b = Vector.new(...)
u, v, s = SV.decomp(m)
x = SV.solve(u, b, s, b)

ex2:
```
x = SV.solve(m, b)
```
ex3:
```
x = m.SV_solve(b)
```

Householder solver for linear systems

GSL::Linalg::HH.solve(A, b)
GSL::Matrix#HH_solve(b): These methods solve the system A x = b directly using Householder transformations.
GSL::Linalg::HH.svx(A, b)
GSL::Matrix#HH_svx(b): These methods solve the system A x = b in-place directly using Householder transformations. The input vector b is replaced by the solution.

back