Basic Linear Algebra

Methods

================================ ============ Method Description ================================ ============ inv_ (A) Inverse solve_ (A, b) Solve equation for general matrix solveBanded_ (kl, ku, ab, b) Solve equation for banded matrix solveHBanded_ (AB, b[, lower]) Solve equation for symmetric banded matrix solveTriangular_ (A, b, lower) Solve equation for triangular matrix det_ (A) Determinant norm_ (A, norm) Norm of a matrix lstsq_ (A, b[, cond]) Least-squares solution pinv_ (A) Moore-Penrose pseudoinverse kron_ (A, B) Kronecker product orthogonalProcrustes_ (A, B) Solution to the orthogonal Procrustes problem subspaceAngles_ (A, B) Subspace angles bandwidth_ (A) Lower and upper bandwidth issymmetric_ (A[, atol]) Check symmetric condition ishermitian_ (A[, atol]) Check hermitian condition ================================ ============

Description

.. _inv:

.. function:: inv(a)

:param a: Matrix<Float64>
:return: Matrix<Float64>

Compute the inverse of a matrix.

.. _solve:

.. function:: solve(A, b)

:param A: Matrix<Float64>
:param b: Matrix<Float64>
:return: Matrix<Float64>

Solve the equation :math:`Ax = b`, where :math:`A` is a general matrix.

.. _solveBanded:

.. function:: solveBanded(kl, ku, A, b)

:param kl: Int64, Number of sub-diagonals in the lower-triangle
:param ku: Int64, Number of sub-diagonals in the upper-triangle
:param A: Matrix<Float64>
:param b: Matrix<Float64>
:return: Matrix<Float64>

Solve the equation :math:`Ax = b`, where :math:`A` is a banded matrix.

.. _solveHBanded:

.. function:: solveHBanded(AB, b[, lower])

:param AB: Matrix<Float64>
:param b: Matrix<Float64>
:param lower: Bool
:return: Matrix<Float64>

Solve the equation :math:`Ax = b`, where :math:`A` is a symmetric banded matrix.

.. _solveTriangular:

.. function:: solveTriangular(A, b, lower)

:param A: Matrix<Float64>
:param b: Matrix<Float64>
:param lower: Bool
:return: Matrix<Float64>

Solve the equation :math:`Ax = b`, where :math:`A` is a triangular matrix.

.. _det:

.. function:: det(A)

:param A: Matrix<Float64>
:return: Float64

Compute the determinant of a matrix.

.. _norm:

.. function:: norm(A, norm)

:param A: Matrix<Float64>
:param norm: Rune, type of norm to compute. One of `F`, `I`, `1` and `M`.
:return: Float64

Compute the norm of a matrix.

.. _lstsq:

.. function:: lstsq(A, b[, cond])

:param A: Matrix<Float64>
:param b: Vector<Float64>
:param cond: Float64. Cutoff for singular values. Optional, default to 0.0. 
:return: Vector<Float64>

Compute least-squares solution to the equation :math:`Ax=b`.

Let :math:`A` be the Moore-Penrose pseudoinverse of :math:`A` (computed by
:ref:`pinv <pinv>`), then the least-squares is :math:`x=A^+b`.

Example:

1. Suppose we wish to fit the following two vectors:

.. math::

    x = [1.0, 2.5, 3.5, 4.0, 5.0, 7.0, 8.5]

    y = [0.3, 1.1, 1.5, 2.0, 3.2, 6.6, 8.6]

to the equation :math:`y=a+bx^2`. Construct the 7-by-2 matrix :math:`M`, whose
first column are 1's, and second column are squares of :math:`x`. Then the
least-squares solution to the equation :math:`Mp=y` gives the coefficients as
vector :math:`p`.

::

    let x = vector([1.0, 2.5, 3.5, 4.0, 5.0, 7.0, 8.5])
    let y = vector([0.3, 1.1, 1.5, 2.0, 3.2, 6.6, 8.6])
    let M = empty<Float64>(x.size(), 2)
    for (i in 0..x.size()) {
        M[i,0] = 1.0
        M[i,1] = pow(x[i], 2.0)
    }
    let p = lstsq(M, y)

Result: :math:`[0.209258, 0.120139]`

.. _pinv:

.. function:: pinv(A[, atol, rtol])

:param A: Matrix<Float64>
:param atol: absolute tolerance
:param rtol: relative tolerance
:return: Matrix<Float64>, the Moore-Penrose pseudoinverse of :math:`A`.

Compute the Moore-Penrose pseudoinverse of :math:`A`.

First, compute the singular value decomposition :math:`A=USV^T`, then
take the :math:`T` to be the diagonal matrix that inverts every nonzero
value on the diagonal of :math:`S`. Here a singular value is determined
to be nonzero if it is greater than ``atol + max(s) * rtol``. If ``rtol``
is set to zero, it is considered to be floating point eps multiplied by
the maximum dimension of :math:`A`.

Example:

1. The pseudoinverse of the following matrix :math:`A` is itself:

.. math::

    \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{bmatrix}

::

    var A = matrix<Float64>(
        [[1.0, 0.0, 0.0],
         [0.0, 1.0, 0.0],
         [0.0, 0.0, 0.0]]
    )
    let pinvA = pinv(A)

.. _kron:

.. function:: kron(A, B)

:param T: type parameter, implements ``Number``
:param A: Matrix<T>
:param B: Matrix<T>
:return: Matrix<T>

Compute the Kronecker product of two matrices.

.. _orthogonalProcrustes:

.. function:: orthogonalProcrustes(A, B)

:param A: Matrix<Float64>
:param B: Matrix<Float64>
:return: (Matrix<Float64>, Float64)

Compute the solution to the orthogonal Procrustes problem.

Find the orthogonal matrix :math:`R` that minimizes the Frobenius norm
of :math:`AR-B`. Let :math:`A^TB=USV^T` be the singular value decomposition
of :math:`A^TB`. Then :math:`R=UV^T`. The returned scale is the sum of
singular values in :math:`S`.

Example:

1. The following problem has an exact solution. Consider two matrices as follows:

.. math::

    A = \begin{bmatrix} 2 & 0 & 1 \\ -2 & 0 & 0 \end{bmatrix}

    B = \begin{bmatrix} 1 & 0 & 2 \\ 0 & 0 & -2 \end{bmatrix}

Then :math:`A` can be transformed to :math:`B` by swapping the first and third
coordinates. This can be verified as follows:

::

    var a = matrix<Float64>(
        [[ 2.0, 0.0, 1.0],
        [-2.0, 0.0, 0.0]]
    )
    var b = matrix<Float64>(
        [[1.0, 0.0,  2.0],
        [0.0, 0.0, -2.0]]
    )
    var (R, sca) = orthogonalProcrustes(a, b)

Result: sca = 9.0 and

.. math::

    R = \begin{bmatrix} 0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{bmatrix}

.. _subspaceAngles:

.. function:: subspaceAngles(A, B)

:param A: Matrix<Float64>
:param B: Matrix<Float64>
:return: Vector<Float64>, the list of angles between subspaces.

Compute the subspace angles between :math:`A` and :math:`B`. The two
input matrices should have the same number of rows.

Example:

1. The angles between first two columns and last two columns of size-4
Hadamard matrix should be :math:`\frac{\pi}{2}`.

::

    var A = matrix<Float64>(
        [[1.0,  1.0],
         [1.0, -1.0],
         [1.0,  1.0],
         [1.0, -1.0]]
    )
    var B = matrix<Float64>(
        [[ 1.0,  1.0],
         [ 1.0, -1.0],
         [-1.0, -1.0],
         [-1.0,  1.0]]
    )
    let angles = subspaceAngles(A, B)

Result: :math:`[1.570796, 1.570796]`

2. The angles between first two columns of size-4 Hadamard matrix and itself
should be zero.

::

    var A = matrix<Float64>(
        [[1.0,  1.0],
         [1.0, -1.0],
         [1.0,  1.0],
         [1.0, -1.0]]
    )
    let angles = subspaceAngles(A, A)

Result: :math:`[0.0, 0.0]`

.. _bandwidth:

.. function:: bandwidth(A)

:param A: Matrix<Float64>
:return: (Int64, Int64), the lower and upper bandwidth.

Compute the lower and upper bandwidth of :math:`A`. This is the number of
non-zero subdiagonals below and above the main diagonal in :math:`A`.

.. _issymmetric:

.. function:: issymmetric(A[, atol])

:param A: Matrix<Float64>
:param atol: Float64, absolute tolerance. Optional, default to 0.
:return: Bool

Return whether the matrix :math:`A` is symmetric up to the given tolerance.

.. _ishermitian:

.. function:: ishermitian(A[, atol])

:param A: Matrix<Complex64>
:param atol: Float64, absolute tolerance. Optional, default to 0.
:return: Bool

Return whether the matrix :math:`A` is hermitian up to the given tolerance.