# Linalg.Generic¶

This document is auto-generated for Owl’s APIs. #66 entries have been extracted.

Github: {Signature} {Implementation}

## 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.)

source code

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.

source code

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


det x computes the determinant of a square matrix x.

source code

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.

source code

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.

source code

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.

source code

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 Dense.Ndarray.Generic, this function assumes the input is either 1d vector or 2d matrix.

$||v||_p = \Big[ \sum_{k=0}^{N-1} |v_k|^p \Big]^{1/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 $$||v||_{\infty} = \max_i(|v(i)|)$$.
• If p = -infinity, then returns $$||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.

source code

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.

source code

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.

source code

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


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

source code

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


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

source code

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


is_symmetric x returns true if x is symmetric otherwise false.

source code

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


is_hermitian x returns true if x is hermitian otherwise false.

source code

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


is_diag x returns true if x is diagonal otherwise false.

source code

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


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

source code

## Factorisation¶

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 lu : ('a, 'b) t -> ('a, 'b) t * ('a, 'b) t * (int32, int32_elt) t


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

source code

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.

source code

val qr : ?thin:bool -> ?pivot:bool -> ('a, 'b) t -> ('a, 'b) t * ('a, 'b) t * (int32, 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.

source code

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.

source code

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.

source code

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.

source code

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.

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

source code

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.

source code

val schur : otyp:('c, 'd) 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.

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

source code

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


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

source code

val ordschur : otyp:('c, 'd) kind -> select:(int32, 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

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

source code

val qz : otyp:('c, 'd) 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.

$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).

source code

val ordqz : otyp:('c, 'd) kind -> select:(int32, 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).

source code

val qzvals : otyp:('c, 'd) 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.

source code

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

$X = Q H Q^T$

source code

## Eigenvalues & eigenvectors¶

val eig : ?permute:bool -> ?scale:bool -> otyp:('a, 'b) 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.

source code

val eigvals : ?permute:bool -> ?scale:bool -> otyp:('a, 'b) 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.

source code

## 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,

$X^T X = I$

source code

val linsolve : ?trans:bool -> ('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. The function uses 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.

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

$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 Operator Functor.

source code

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


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

$Y = A + BX$

source code

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.

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

source code

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.

$AX + XA^H = Q$
Parameters:
• a : m x m matrix A.
• q : n x n matrix Q.
Returns:
• x : m x n matrix X.

source code

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

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

source code

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

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

source code

## Low-level factorisation functions¶

val lufact : ('a, 'b) t -> ('a, 'b) t * (int32, int32_elt) t


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

source code

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


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

source code

val bkfact : ?upper:bool -> ?symmetric:bool -> ?rook:bool -> ('a, 'b) t -> ('a, 'b) t * (int32, 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.

source code

## 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.)

source code

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


expm x computes the matrix exponential of x defined by

$e^x = \sum_{k=0}^{\infty} \frac{1}{k!} x^k$

The function implements the scaling and squaring algorithm which uses Padé approximation to compute the matrix exponential [AMH09].

source code

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.

source code

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.

source code

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.

source code

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


sincosm x returns both matrix sine and cosine of x.

source code

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.

source code

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.

source code

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.

source code

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


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

source code

## Helper functions¶

val select_ev : [ LHP | RHP | UDI | UDO ] -> ('a, 'b) t -> (int32, 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 $$(real(e) < 0)$$.
• RHP: Left-half plane $$(real(e) \ge 0)$$.
• UDI: Left-half plane $$(abs(e) < 1)$$.
• UDO: Left-half plane $$(abs(e) \ge 0)$$.

source code

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.

source code