Title: Linear Algebra

0.1 Contents:

  1. LU Decomposition
  2. Complex LU Decomposition
  3. QR Decomposition
  4. QR Decomposition with Column Pivoting
  5. Singular Value Decomposition
  6. holesky Decomposition
  7. Bidiagonalization
  8. Householder Transformations
  9. Householder solver for linear systems
  10. Tridiagonal Systems
  11. Balancing
  12. NArray

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

Examples:

  1. Singleton method of the GSL::Linalg::LU module

    irb(main):012:0> m = Matrix[1..9, 3, 3]
    => GSL::Matrix: 
    [ 1.000e+00 2.000e+00 3.000e+00 
      4.000e+00 5.000e+00 6.000e+00 
      7.000e+00 8.000e+00 9.000e+00 ]
    irb(main):013:0> lu, perm, sign = Linalg::LU.decomp(m)
  2. Instance method of GSL::Matrix class

    irb(main):012:0> 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::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.

  1. Singleton method of the GSL::Linalg::LU module

    A = Matrix[[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         
    b = Vector[1, 2, 3, 4]
    x = Linalg::LU.solve(lu, perm, b)    
  2. Instance method of GSL::Linalg::LUMatrix class

    lu, perm, sign = A.LU_decomp  # lu is an instance of Linalg::LUMatrix class
    b = Vector[1, 2, 3, 4]
    x = lu.solve(perm, b)    
  3. Solve directly

    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::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::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::LUMatrix#det(signum)
These methods return the determinant of the matrix.

2 Complex LU decomposition

3 QR decomposition

GSL::Linalg::QR_decomp(A)
GSL::Linalg::QR.decomp(A)
GSL::Matrix#QR_decomp
These compute QR decomposition of the matrix and return an array [QR, tau].
  1. Singleton method of the module GSL::Linalg

    qr, tau = Linalg::QR_decomp(m)
    p qr.class                 # GSL::Linalg::QRMatrix, subclass of GSL::Matrix
    p tau.class                # GSL::Linalg::TauVector, subclass of GSL::Vector
  2. Singleton method of the module GSL::Linalg:QR

    qr, tau = Linalg::QR.decomp(m)
  3. Instance method of GSL::Matrix

    qr, tau = m.QR_decomp
GSL::Linalg::QR.solve(A, b)
GSL::Linalg::QR.solve(QR, tau, b)
GSL::Matrix#QR_solve(b)
GSL::Linalg::QRMatrix#solve(tau, b)
Solve the system A x = b using the QR decomposition.
GSL::Linalg::QR.svx(A, x)
GSL::Linalg::QR.svx(QR, tau, x)
GSL::Matrix#QR_svx(x)
GSL::Linalg::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::QRMatrix#unpack(tau)

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

Ex:

irb(main):010:0> m = Matrix[1..9, 3, 3]
=> GSL::Matrix: 
[ 1.000e+00 2.000e+00 3.000e+00 
  4.000e+00 5.000e+00 6.000e+00 
  7.000e+00 8.000e+00 9.000e+00 ]
irb(main):011:0> qr, tau = m.QR_decomp
irb(main):012:0> q, r = qr.unpack(tau)  
irb(main):013:0> q*r               # Reconstruct the metrix m
=> GSL::Matrix: 
[ 1.000e+00 2.000e+00 3.000e+00 
  4.000e+00 5.000e+00 6.000e+00 
  7.000e+00 8.000e+00 9.000e+00 ]
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.

4 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].
GSL::Linalg::QRPT.decomp2(A)
GSL::Matrix#QRPT_decomp2
These return an array [Q, R, tau, perm, signum].
GSL::Linalg::QRPT.solve(m, b)
GSL::Linalg::QRPT.solve(qr, tau, perm, b)
GSL::Matrix#QRPT_solve(A, b)
GSL::Linalg::QRPQMatrix#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.
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.
GSL::Linalg::QRPT.update(q, r, perm, u, v)
GSL::Linalg::QRPT.Rsolve(qr, perm, b)
GSL::Linalg::QRPTMatrix#Rsolve(perm, b)
GSL::Linalg::QRPT.Rsvx(qr, perm, b)
GSL::Linalg::QRPTMatrix#Rsvx(perm, b)

5 Singular Value Decomposition

GSL::Linalg::SV.decomp(A[, work])
GSL::Matrix#SV_decomp([work])

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

Ex:

irb(main):020:0* m = Matrix[[3, 5, 2], [5, 1, 4], [7, 6, 3]]
=> GSL::Matrix: 
[ 3.000e+00 5.000e+00 2.000e+00 
  5.000e+00 1.000e+00 4.000e+00 
  7.000e+00 6.000e+00 3.000e+00 ]
irb(main):021:0> u, v, s = m.SV_decomp  # u, v: Matrix, s: Vector (singular values)
irb(main):022:0> u*u.trans              # u is orthnormal
=> GSL::Matrix: 
[  1.000e+00  2.452e-17 -4.083e-16 
   2.452e-17  1.000e+00 -3.245e-16 
  -4.083e-16 -3.245e-16  1.000e+00 ]
irb(main):023:0> v*v.trans              # v is also orthnormal
=> GSL::Matrix: 
[  1.000e+00  3.555e-17 -1.867e-16 
   3.555e-17  1.000e+00 -1.403e-16 
  -1.867e-16 -1.403e-16  1.000e+00 ]
irb(main):024:0> u*Matrix.diagonal(s)*v.trans # Reconstruct the matrix
=> GSL::Matrix: 
[ 3.000e+00 5.000e+00 2.000e+00 
  5.000e+00 1.000e+00 4.000e+00 
  7.000e+00 6.000e+00 3.000e+00 ]
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.

6 Cholesky Decomposition

A symmetric, positive definite square matrix A has a Cholesky decomposition into a product of a lower triangular matrix L and its transpose L^T, as A = L L^T. This is sometimes referred to as taking the square-root of a matrix. The Cholesky decomposition can only be carried out when all the eigenvalues of the matrix are positive. This decomposition can be used to convert the linear system A x = b into a pair of triangular systems (L y = b, L^T x = y), which can be solved by forward and back-substitution.

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

This method factorizes the positive-definite square matrix A into the Cholesky decomposition A = L L^T. The upper triangular part of the matrix returned contains L^T, the diagonal terms being identical for both L and L^T. If the matrix is not positive-definite then the decomposition will fail.

Ex:

irb(main):006:0> m = Matrix.pascal(3)
=> GSL::Matrix
[  1.000e+00  1.000e+00  1.000e+00 
   1.000e+00  2.000e+00  3.000e+00 
   1.000e+00  3.000e+00  6.000e+00 ]
irb(main):007:0> c = Linalg::Cholesky.decomp(m)
=> GSL::Linalg::Cholesky::CholeskyMatrix
[  1.000e+00  1.000e+00  1.000e+00 
   1.000e+00  1.000e+00  2.000e+00 
   1.000e+00  2.000e+00  1.000e+00 ]
irb(main):008:0> l = c.lower
=> GSL::Matrix
[  1.000e+00  0.000e+00  0.000e+00 
   1.000e+00  1.000e+00  0.000e+00 
   1.000e+00  2.000e+00  1.000e+00 ]
irb(main):009:0> (l*l.trans) == m
=> true
GSL::Linalg::Cholesky.solve(cholesky, b)
GSL::Linalg::Cholesky.svx(cholesky, x)
These methods solve the system A x = b using the Cholesky decomposition of A into the matrix cholesky given by GSL::Linalg::Cholesky.decomp.

7 Bidiagonalization

GSL::Linalg::Bidiag::decomp!(A)
GSL::Linalg::bidiag_decomp!(A)
GSL::Linalg::Bidiag::decomp(A)
GSL::Linalg::bidiag_decomp(A)
GSL::Linalg::Bidiag::unpack
GSL::Linalg::bidiag_unpack
GSL::Linalg::Bidiag::unpack2
GSL::Linalg::bidiag_unpack2
GSL::Linalg::Bidiag::unpack_B
GSL::Linalg::bidiag_unpack_B

8 Householder Transformations

GSL::Linalg::Householder::transform(v)
GSL::Linalg::HH::transform(v)
GSL::Vector#householder_transform
These methods prepare a Householder transformation P = I - tau v v^T which can be used to zero all the elements of the input vector except the first. On output the transformation is stored in the vector v and the scalar tau is returned.
GSL::Linalg::Householder::hm(tau, v, A)
GSL::Linalg::HH::hm(tau, v, A)
These methods apply the Householder matrix P defined by the scalar tau and the vector v to the left-hand side of the matrix A. On output the result P A is stored in A.
GSL::Linalg::Householder::mh(tau, v, A)
GSL::Linalg::HH::mh(tau, v, A)
These methods apply the Householder matrix P defined by the scalar tau and the vector v to the right-hand side of the matrix A. On output the result A P is stored in A.
GSL::Linalg::Householder::hv(tau, v, w)
GSL::Linalg::HH::hv(tau, v, w)
These methods apply the Householder transformation P defined by the scalar tau and the vector v to the vector w. On output the result P w is stored in w.

9 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. The matrix A is not modified.
GSL::Linalg::HH.solve!(A, b)
GSL::Matrix#HH_solve!(b)
These methods solve the system A x = b directly using Householder transformations. The matrix A is destroyed by the 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.

10 Tridiagonal Systems

GSL::Linglg::solve_tridiag(diag, e, f, b)
GSL::Linglg::Tridiag::solve(diag, e, f, b)

These methods solve the general N-by-N system A x = b where A is tridiagonal ( N >= 2). The super-diagonal and sub-diagonal vectors e and f must be one element shorter than the diagonal vector diag. The form of A for the 4-by-4 case is shown below,

A = ( d_0 e_0  0   0  )
    ( f_0 d_1 e_1  0  )
    (  0  f_1 d_2 e_2 )
    (  0   0  f_2 d_3 )
GSL::Linglg::solve_symm_tridiag(diag, e, b)
GSL::Linglg::Tridiag::solve_symm(diag, e, b)

These methods solve the general N-by-N system A x = b where A is symmetric tridiagonal ( N >= 2). The off-diagonal vector e must be one element shorter than the diagonal vector diag. The form of A for the 4-by-4 case is shown below,

A = ( d_0 e_0  0   0  )
    ( e_0 d_1 e_1  0  )
    (  0  e_1 d_2 e_2 )
    (  0   0  e_2 d_3 )
GSL::Linglg::solve_cyc_tridiag(diag, e, f, b)
GSL::Linglg::Tridiag::solve_cyc(diag, e, f, b)

These methods solve the general N-by-N system A x = b where A is cyclic tridiagonal ( N >= 3). The cyclic super-diagonal and sub-diagonal vectors e and f must have the same number of elements as the diagonal vector diag. The form of A for the 4-by-4 case is shown below,

A = ( d_0 e_0  0  f_3 )
    ( f_0 d_1 e_1  0  )
    (  0  f_1 d_2 e_2 )
    ( e_3  0  f_2 d_3 )
GSL::Linglg::solve_symm_cyc_tridiag(diag, e, b)
GSL::Linglg::Tridiag::solve_symm_cyc(diag, e, b)

These methods solve the general N-by-N system A x = b where A is symmetric cyclic tridiagonal ( N >= 3). The cyclic off-diagonal vector e must have the same number of elements as the diagonal vector diag. The form of A for the 4-by-4 case is shown below,

A = ( d_0 e_0  0  e_3 )
    ( e_0 d_1 e_1  0  )
    (  0  e_1 d_2 e_2 )
    ( e_3  0  e_2 d_3 )

11 Balancing

GSL::Linalg.balance_columns(A, D)
GSL::Matrix#balance_columns(D)

12 NArray

The following Linalg methods can handle NArray objects:

prev next

Reference index top