# Module Owl_linalg_generic

Linear algebra module including high-level functions to solve linear systems, factorisation, and etc.

###### 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, Stdlib.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, Stdlib.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) Stdlib.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) Stdlib.Bigarray.kind -> select:(int32, Stdlib.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) Stdlib.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) Stdlib.Bigarray.kind -> select:(int32, Stdlib.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) Stdlib.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) Stdlib.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) Stdlib.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 : ?⁠diag_r:bool -> (float, 'a) t -> (float, 'a) t -> (float, 'a) t -> (float, 'a) t -> (float, 'a) t

care ?diag_r 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. * diag_r : true if R is a diagonal matrix. false by default

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, Stdlib.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, Stdlib.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, Stdlib.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 rth 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, Stdlib.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.