package owl

  1. Overview
  2. Docs
package owl
Legend:
Library
Module
Module type
Parameter
Class
Class type
include module type of struct include Owl_linalg_generic end

The module includes a set of advanced linear algebra operations such as singular value decomposition, and etc.

Currently, Linalg module supports dense matrix of four different number types, including ``float32``, ``float64``, ``complex32``, and ``complex64``. The support for sparse matrices will be provided in future.

Type definition
type ('a, 'b) t = ('a, 'b) Owl_dense_matrix_generic.t

Matrix type, a special case of N-dimensional array.

Basic functions
val inv : ('a, 'b) t -> ('a, 'b) t

``inv x`` calculates the inverse of an invertible square matrix ``x`` such that ``x *@ x = I`` wherein ``I`` is an identity matrix. (If ``x`` is singular, ``inv`` will return a useless result.)

val pinv : ?tol:float -> ('a, 'b) t -> ('a, 'b) t

``pinv x`` computes Moore-Penrose pseudoinverse of matrix ``x``. ``tol`` specifies the tolerance, the absolute value of the elements smaller than ``tol`` will be set to zeros.

val det : ('a, 'b) t -> 'a

``det x`` computes the determinant of a square matrix ``x``.

val logdet : ('a, 'b) t -> 'a

``logdet x`` computes the log of the determinant of a square matrix ``x``. It is equivalent to ``log (det x)`` but may provide more accuracy and efficiency.

val rank : ?tol:float -> ('a, 'b) t -> int

``rank x`` calculates the rank of a rectangular matrix ``x`` of shape ``m x n``. The function does so by counting the number of singular values of ``x`` which are beyond a pre-defined threshold ``tol``. By default, ``tol = max(m,n) * eps`` where ``eps = 1e-10``.

val norm : ?p:float -> ('a, 'b) t -> float

``norm ~p x`` computes the matrix p-norm of the passed in matrix ``x``.

Parameters: * ``p`` is the order of norm, the default value is 2. * ``x`` is the input matrix.

Returns: * If ``p = 1``, then returns the maximum absolute column sum of the matrix. * If ``p = 2``, then returns approximately ``max (svd x)``. * If ``p = infinity``, then returns the maximum absolute row sum of the matrix. * If ``p = -1``, then returns the minimum absolute column sum of the matrix. * If ``p = -2``, then returns approximately ``min (svd x)``. * If ``p = -infinity``, then returns the minimum absolute row sum of the matrix.

val vecnorm : ?p:float -> ('a, 'b) t -> float

``vecnorm ~p x`` calculates the generalised vector p-norm, defined as below. If ``x`` is a martrix, it will be flatten to a vector first. Different from the function of the same name in :doc:`owl_dense_ndarray_generic`, this function assumes the input is either 1d vector or 2d matrix.

.. math:: ||v||_p = \Big \sum_{k=0}^{N-1} |v_k|^p \Big^

/p

Parameters: * ``p`` is the order of norm, the default value is 2. * ``x`` is the input vector or matrix.

Returns: * If ``p = infinity``, then returns :math:`||v||_\infty = \max_i(|v(i)|)`. * If ``p = -infinity``, then returns :math:`||v||_

\infty

}

= \min_i(|v(i)|)`. * If ``p = 2`` and ``x`` is a matrix, then returns Frobenius norm of ``x``. * Otherwise returns generalised vector p-norm defined above.

val cond : ?p:float -> ('a, 'b) t -> float

``cond ~p x`` computes the p-norm condition number of matrix ``x``.

``cond ~p:1. x`` returns the 1-norm condition number;

``cond ~p:2. x`` or ``cond x`` returns the 2-norm condition number.

``cond ~p:infinity x`` returns the infinity norm condition number.

The default value of ``p`` is ``2.``

val rcond : ('a, 'b) t -> float

``rcond x`` returns an estimate for the reciprocal condition of ``x`` in 1-norm. If ``x`` is well conditioned, the returned result is near ``1.0``. If ``x`` is badly conditioned, the result is near ``0.``

Check matrix types
val is_square : ('a, 'b) t -> bool

``is_square x`` returns ``true`` if ``x`` is a square matrix otherwise ``false``.

val is_triu : ('a, 'b) t -> bool

``is_triu x`` returns ``true`` if ``x`` is upper triangular otherwise ``false``.

val is_tril : ('a, 'b) t -> bool

``is_tril x`` returns ``true`` if ``x`` is lower triangular otherwise ``false``.

val is_symmetric : ('a, 'b) t -> bool

``is_symmetric x`` returns ``true`` if ``x`` is symmetric otherwise ``false``.

val is_hermitian : (Stdlib.Complex.t, 'a) t -> bool

``is_hermitian x`` returns ``true`` if ``x`` is hermitian otherwise ``false``.

val is_diag : ('a, 'b) t -> bool

``is_diag x`` returns ``true`` if ``x`` is diagonal otherwise ``false``.

val is_posdef : ('a, 'b) t -> bool

``is_posdef x`` checks whether ``x`` is a positive semi-definite matrix.

Factorisation
val lu : ('a, 'b) t -> ('a, 'b) t * ('a, 'b) t * (int32, Bigarray.int32_elt) t

``lu x -> (l, u, ipiv)`` calculates LU decomposition of ``x``. The pivoting is used by default.

val lq : ?thin:bool -> ('a, 'b) t -> ('a, 'b) t * ('a, 'b) t

``lq x -> (l, q)`` calculates the LQ decomposition of ``x``. By default, the reduced LQ decomposition is performed. But you can get full ``Q`` by setting parameter ``thin = false``.

val qr : ?thin:bool -> ?pivot:bool -> ('a, 'b) t -> ('a, 'b) t * ('a, 'b) t * (int32, Bigarray.int32_elt) t

``qr x`` calculates QR decomposition for an ``m`` by ``n`` matrix ``x`` as ``x = Q R``. ``Q`` is an ``m`` by ``n`` matrix (where ``Q^T Q = I``) and ``R`` is an ``n`` by ``n`` upper-triangular matrix.

The function returns a 3-tuple, the first two are ``q`` and ``r``, and the thrid is the permutation vector of columns. The default value of ``pivot`` is ``false``, setting ``pivot = true`` lets ``qr`` performs pivoted factorisation. Note that the returned indices are not adjusted to 0-based C layout.

By default, ``qr`` performs a reduced QR factorisation, full factorisation can be enabled by setting ``thin`` parameter to ``false``.

val chol : ?upper:bool -> ('a, 'b) t -> ('a, 'b) t

``chol x -> u`` calculates the Cholesky factorisation of a positive definite matrix ``x`` such that ``x = u' *@ u``. By default, the upper triangular matrix is returned. The lower triangular part can be obtained by setting the parameter ``upper = false``.

val svd : ?thin:bool -> ('a, 'b) t -> ('a, 'b) t * ('a, 'b) t * ('a, 'b) t

``svd x -> (u, s, vt)`` calculates the singular value decomposition of ``x``, and returns a 3-tuple ``(u,s,vt)``. By default, a reduced svd is performed: E.g., for a ``m x n`` matrix ``x`` wherein ``m <= n``, ``u`` is returned as an ``m`` by ``m`` orthogonal matrix, ``s`` an ``1`` by ``m`` row vector of singular values, and ``vt`` is the transpose of an ``n`` by ``m`` orthogonal rectangular matrix.

The full svd can be performed by setting ``thin = false``. Note that for complex numbers, the type of returned singular values are also complex, the imaginary part is zero.

val svdvals : ('a, 'b) t -> ('a, 'b) t

``svdvals x -> s`` performs the singular value decomposition of ``x`` like ``svd x``, but the function only returns the singular values without ``u`` and ``vt``. Note that for complex numbers, the return is also complex type.

val gsvd : ('a, 'b) t -> ('a, 'b) t -> ('a, 'b) t * ('a, 'b) t * ('a, 'b) t * ('a, 'b) t * ('a, 'b) t * ('a, 'b) t

``gsvd x y -> (u, v, q, d1, d2, r)`` computes the generalised singular value decomposition of a pair of general rectangular matrices ``x`` and ``y``. ``d1`` and ``d2`` contain the generalised singular value pairs of ``x`` and ``y``. The shape of ``x`` is ``m x n`` and the shape of ``y`` is ``p x n``.

.. code-block:: ocaml

let x = Mat.uniform 5 5;; let y = Mat.uniform 2 5;; let u, v, q, d1, d2, r = Linalg.gsvd x y;; Mat.(u *@ d1 *@ r *@ transpose q =~ x);; Mat.(v *@ d2 *@ r *@ transpose q =~ y);;

Please refer to: `Intel MKL Reference <https://software.intel.com/en-us/mkl-developer-reference-c-ggsvd3>`_

val gsvdvals : ('a, 'b) t -> ('a, 'b) t -> ('a, 'b) t

``gsvdvals x y`` is similar to ``gsvd x y`` but only returns the singular values of the generalised singular value decomposition of ``x`` and ``y``.

val schur : otyp:('c, 'd) Bigarray.kind -> ('a, 'b) t -> ('a, 'b) t * ('a, 'b) t * ('c, 'd) t

``schur x -> (t, z, w)`` calculates Schur factorisation of ``x`` in the following form.

.. math:: X = Z T Z^H

Parameters: * ``otyp``: the complex type of eigen values. * ``x``: the ``n x n`` square matrix.

Returns: * ``t`` is (quasi) triangular Schur factor. * ``z`` is orthogonal/unitary Schur vectors. The eigen values are not sorted, they have the same order as that they appear on the diagonal of the output of Schur form ``t``. * ``w`` contains the eigen values of ``x``. ``otyp`` is used to specify the type of ``w``. It needs to be consistent with input type. E.g., if the input ``x`` is ``float32`` then ``otyp`` must be ``complex32``. However, if you use S, D, C, Z module, then you do not need to worry about ``otyp``.

val schur_tz : ('a, 'b) t -> ('a, 'b) t * ('a, 'b) t

``schur_tz x`` is similar to ``schur`` but only returns ``(t, z)``.

val ordschur : otyp:('c, 'd) Bigarray.kind -> select:(int32, Bigarray.int32_elt) t -> ('a, 'b) t -> ('a, 'b) t -> ('a, 'b) t * ('a, 'b) t * ('c, 'd) t

``ordschur ~select t z -> (r, p)`` reorders ``t`` and ``z`` returned by Schur factorization ``schur x -> (t, z)`` according ``select`` such that

.. math:: X = P R P^H

Parameters: * ``otyp``: the complex type of eigen values * ``select`` the logical vector to select eigenvalues, refer to ``select_ev``. * ``t``: the Schur matrix returned by ``schur x``. * ``z``: the unitary matrix ``z`` returned by ``schur x``.

Returns: * ``r``: reordered Schur matrix ``t``. * ``p``: reordered orthogonal matrix ``z``.

val qz : otyp:('c, 'd) Bigarray.kind -> ('a, 'b) t -> ('a, 'b) t -> ('a, 'b) t * ('a, 'b) t * ('a, 'b) t * ('a, 'b) t * ('c, 'd) t

``qz x -> (s, t, q, z, w)`` calculates generalised Schur factorisation of ``x`` in the following form. It is also known as QZ decomposition.

.. math:: X = Q S Z^H Y = Z T Z^H

Parameters: * ``otyp``: the complex type of eigen values. * ``x``: the ``n x n`` square matrix. * ``y``: the ``n x n`` square matrix.

Returns: * ``s``: the upper quasitriangular matrices S. * ``t``: the upper quasitriangular matrices T. * ``q``: the unitary matrices Q. * ``z``: the unitary matrices Z. * ``w``: the generalised eigenvalue for a pair of matrices (X,Y).

val ordqz : otyp:('c, 'd) Bigarray.kind -> select:(int32, Bigarray.int32_elt) t -> ('a, 'b) t -> ('a, 'b) t -> ('a, 'b) t -> ('a, 'b) t -> ('a, 'b) t * ('a, 'b) t * ('a, 'b) t * ('a, 'b) t * ('c, 'd) t

``ordqz ~select a b q z`` reorders the generalised Schur decomposition of a pair of matrices (X,Y) so that a selected cluster of eigenvalues appears in the leading diagonal blocks of (X,Y).

val qzvals : otyp:('c, 'd) Bigarray.kind -> ('a, 'b) t -> ('a, 'b) t -> ('c, 'd) t

``qzvals ~otyp x y`` is similar to ``qz ~otyp x y`` but only returns the generalised eigen values.

val hess : ('a, 'b) t -> ('a, 'b) t * ('a, 'b) t

``hess x -> (h, q)`` calculates the Hessenberg form of a given matrix ``x``. Both Hessenberg matrix ``h`` and unitary matrix ``q`` is returned, such that ``x = q *@ h *@ (transpose q)``.

.. math:: X = Q H Q^T

Eigenvalues & eigenvectors
val eig : ?permute:bool -> ?scale:bool -> otyp:('a, 'b) Bigarray.kind -> ('c, 'd) t -> ('a, 'b) t * ('a, 'b) t

``eig x -> v, w`` computes the right eigenvectors ``v`` and eigenvalues ``w`` of an arbitrary square matrix ``x``. The eigenvectors are column vectors in ``v``, their corresponding eigenvalues have the same order in ``w`` as that in ``v``.

Note that ``otyp`` specifies the complex type of the output, but you do not need worry about this parameter if you use S, D, C, Z modules in Linalg.

val eigvals : ?permute:bool -> ?scale:bool -> otyp:('a, 'b) Bigarray.kind -> ('c, 'd) t -> ('a, 'b) t

``eigvals x -> w`` is similar to ``eig`` but only computes the eigenvalues of an arbitrary square matrix ``x``.

Linear system of equations
val null : ('a, 'b) t -> ('a, 'b) t

``null a -> x`` computes an orthonormal basis ``x`` for the null space of ``a`` obtained from the singular value decomposition. Namely, ``a *@ x`` has negligible elements, ``M.col_num x`` is the nullity of ``a``, and ``transpose x *@ x = I``. Namely,

.. math:: X^T X = I

val triangular_solve : upper:bool -> ?trans:bool -> ('a, 'b) t -> ('a, 'b) t -> ('a, 'b) t

``triangular_linsolve a b -> x`` solves a linear system of equations ``a * x = b`` where ``a`` is either a upper or a lower triangular matrix. This function uses cblas ``trsm`` under the hood.

.. math:: AX = B

By default, ``trans = false`` indicates no transpose. If ``trans = true``, then function will solve ``A^T * x = b`` for real matrices; ``A^H * x = b`` for complex matrices.

.. math:: A^H X = B

val linsolve : ?trans:bool -> ?typ:[ `n | `u | `l ] -> ('a, 'b) t -> ('a, 'b) t -> ('a, 'b) t

``linsolve a b -> x`` solves a linear system of equations ``a * x = b`` in the following form. By default, ``typ=`n`` and the function use LU factorisation with partial pivoting when ``a`` is square and QR factorisation with column pivoting otherwise. The number of rows of ``a`` must equal the number of rows of ``b``. If ``a`` is a upper(lower) triangular matrix, the function calls the ``solve_triangular`` function when ``typ=`u``(``typ=`l``).

.. math:: AX = B

By default, ``trans = false`` indicates no transpose. If ``trans = true``, then function will solve ``A^T * x = b`` for real matrices; ``A^H * x = b`` for complex matrices.

.. math:: A^H X = B

The associated operator is ``/@``, so you can simply use ``a /@ b`` to solve the linear equation system to get ``x``. Please refer to :doc:`owl_operator`.

val linreg : ('a, 'b) t -> ('a, 'b) t -> 'a * 'a

``linreg x y -> (a, b)`` solves ``y = a + b*x`` using Ordinary Least Squares.

.. math:: Y = A + BX

val sylvester : ('a, 'b) t -> ('a, 'b) t -> ('a, 'b) t -> ('a, 'b) t

``sylvester a b c`` solves a Sylvester equation in the following form. The function calls LAPACKE function ``trsyl`` solve the system.

.. math:: AX + XB = C

Parameters: * ``a`` : ``m x m`` matrix A. * ``b`` : ``n x n`` matrix B. * ``c`` : ``m x n`` matrix C.

Returns: * ``x`` : ``m x n`` matrix X.

val lyapunov : ('a, 'b) t -> ('a, 'b) t -> ('a, 'b) t

``lyapunov a q`` solves a continuous Lyapunov equation in the following form. The function calls LAPACKE function ``trsyl`` solve the system. In Matlab, the same function is called ``lyap``.

.. math:: AX + XA^H = Q

Parameters: * ``a`` : ``m x m`` matrix A. * ``q`` : ``n x n`` matrix Q.

Returns: * ``x`` : ``m x n`` matrix X.

val discrete_lyapunov : ?solver:[ `default | `bilinear | `direct ] -> ('a, 'b) t -> ('a, 'b) t -> ('a, 'b) t

``discrete_lyapunov a q`` solves a discrete-time Lyapunov equation in the following form.

.. math:: X - AXA^H = Q

Parameters: * ``a`` : ``m x m`` matrix A. * ``q`` : ``n x n`` matrix Q.

Returns: * ``x`` : ``m x n`` matrix X.

val care : (float, 'a) t -> (float, 'a) t -> (float, 'a) t -> (float, 'a) t -> (float, 'a) t

``care a b q r`` solves the continuous-time algebraic Riccati equation system in the following form. The algorithm is based on :cite:`laub1979schur`.

.. math:: A^T X + X A − X B R^

1

}

B^T X + Q = 0

Parameters: * ``a`` : real cofficient matrix A. * ``b`` : real cofficient matrix B. * ``q`` : real cofficient matrix Q. * ``r`` : real cofficient matrix R. R must be non-singular.

Returns: * ``x`` : a solution matrix X.

val dare : (float, 'a) t -> (float, 'a) t -> (float, 'a) t -> (float, 'a) t -> (float, 'a) t

``dare a b q r`` solves the discrete-time algebraic Riccati equation system in the following form. The algorithm is based on :cite:`laub1979schur`.

.. math:: A^T X A - X - (A^T X B) (B^T X B + R)^

1

}

(B^T X A) + Q = 0

Parameters: * ``a`` : real cofficient matrix A. A must be non-singular. * ``b`` : real cofficient matrix B. * ``q`` : real cofficient matrix Q. * ``r`` : real cofficient matrix R. R must be non-singular.

Returns: * ``x`` : a symmetric solution matrix X.

Low-level factorisation functions
val lufact : ('a, 'b) t -> ('a, 'b) t * (int32, Bigarray.int32_elt) t

``lufact x -> (a, ipiv)`` calculates LU factorisation with pivot of a general matrix ``x``.

val qrfact : ?pivot:bool -> ('a, 'b) t -> ('a, 'b) t * ('a, 'b) t * (int32, Bigarray.int32_elt) t

``qrfact x -> (a, tau, jpvt)`` calculates QR factorisation of a general matrix ``x``.

val bkfact : ?upper:bool -> ?symmetric:bool -> ?rook:bool -> ('a, 'b) t -> ('a, 'b) t * (int32, Bigarray.int32_elt) t

``bk x -> (a, ipiv)`` calculates Bunch-Kaufman factorisation of ``x``. If ``symmetric = true`` then ``x`` is symmetric, if ``symmetric = false`` then ``x`` is hermitian. If ``rook = true`` the function performs bounded Bunch-Kaufman ("rook") diagonal pivoting method, if ``rook = false`` then Bunch-Kaufman diagonal pivoting method is used. ``a`` contains details of the block-diagonal matrix ``d`` and the multipliers used to obtain the factor ``u`` (or ``l``).

The ``upper`` indicates whether the upper or lower triangular part of ``x`` is stored and how ``x`` is factored. If ``upper = true`` then upper triangular part is stored: ``x = u*d*u'`` else ``x = l*d*l'``.

For ``ipiv``, it indicates the details of the interchanges and the block structure of ``d``. Please refer to the function ``sytrf``, ``hetrf`` in MKL documentation for more details.

Matrix functions
val mpow : ('a, 'b) t -> float -> ('a, 'b) t

``mpow x r`` returns the dot product of square matrix ``x`` with itself ``r`` times, and more generally raises the matrix to the ``r``th power. ``r`` is a float that must be equal to an integer; it can be be negative, zero, or positive. Non-integer exponents are not yet implemented. (If ``r`` is negative, ``mpow`` calls ``inv``, and warnings in documentation for ``inv`` apply.)

val expm : ('a, 'b) t -> ('a, 'b) t

``expm x`` computes the matrix exponential of ``x`` defined by

.. math:: e^x = \sum_k=0^\infty \frac

k! x^k

The function implements the scaling and squaring algorithm which uses Padé approximation to compute the matrix exponential :cite:`al2009new`.

val sinm : ('a, 'b) t -> ('a, 'b) t

``sinm x`` computes the matrix sine of input ``x``. The function uses ``expm`` to compute the matrix exponentials.

val cosm : ('a, 'b) t -> ('a, 'b) t

``cosm x`` computes the matrix cosine of input ``x``. The function uses ``expm`` to compute the matrix exponentials.

val tanm : ('a, 'b) t -> ('a, 'b) t

``tanm x`` computes the matrix tangent of input ``x``. The function uses ``expm`` to compute the matrix exponentials.

val sincosm : ('a, 'b) t -> ('a, 'b) t * ('a, 'b) t

``sincosm x`` returns both matrix sine and cosine of ``x``.

val sinhm : ('a, 'b) t -> ('a, 'b) t

``sinhm x`` computes the hyperbolic matrix sine of input ``x``. The function uses ``expm`` to compute the matrix exponentials.

val coshm : ('a, 'b) t -> ('a, 'b) t

``coshm x`` computes the hyperbolic matrix cosine of input ``x``. The function uses ``expm`` to compute the matrix exponentials.

val tanhm : ('a, 'b) t -> ('a, 'b) t

``tanhm x`` computes the hyperbolic matrix tangent of input ``x``. The function uses ``expm`` to compute the matrix exponentials.

val sinhcoshm : ('a, 'b) t -> ('a, 'b) t * ('a, 'b) t

``sinhcoshm x`` returns both hyperbolic matrix sine and cosine of ``x``.

Helper functions
val select_ev : [ `LHP | `RHP | `UDI | `UDO ] -> ('a, 'b) t -> (int32, Bigarray.int32_elt) t

``select_ev keyword ev`` generates a logical vector (of same shape as ``ev``) from eigen values ``ev`` according to the passed in keywards.

  • ``LHP``: Left-half plane :math:`(real(e) < 0)`.
  • ``RHP``: Left-half plane :math:`(real(e) \ge 0)`.
  • ``UDI``: Left-half plane :math:`(abs(e) < 1)`.
  • ``UDO``: Left-half plane :math:`(abs(e) \ge 0)`.
val peakflops : ?n:int -> unit -> float

``peakflops ()`` returns the peak number of float point operations using ``Owl_cblas_basic.dgemm`` function. The default matrix size is ``2000 x 2000``, but you can change this by setting ``n`` to other numbers as you like.

OCaml

Innovation. Community. Security.

GitHub Discord Twitter Peertube RSS
About
Changelog Releases Industrial Users Academic Users Why OCaml
Ecosystem
Packages Community Events OCaml Planet Jobs
Resources
Install OCaml Get Started Platform Tools Language Manual Standard Library API Books Exercises Papers OCaml Playground Logo
Policies
Carbon Footprint Governance Privacy Code of Conduct