arb_mat.h – matrices over the real numbers¶
An arb_mat_t represents a dense matrix over the real numbers,
implemented as an array of entries of type arb_struct.
The dimension (number of rows and columns) of a matrix is fixed at
initialization, and the user must ensure that inputs and outputs to
an operation have compatible dimensions. The number of rows or columns
in a matrix can be zero.
Note
Methods prefixed with arb_mat_approx treat all input entries as floating-point numbers (ignoring the radii of the balls) and compute floating-point output (balls with zero radius) representing approximate solutions without error bounds. All other methods compute rigorous error bounds. The approx methods are typically useful for computing initial values or preconditioners for rigorous solvers. Some users may also find approx methods useful for doing ordinary numerical linear algebra in applications where error bounds are not needed.
Types, macros and constants¶
-
type arb_mat_struct¶
-
type arb_mat_t¶
Contains a pointer to a flat array of the entries (entries), an array of pointers to the start of each row (rows), and the number of rows (r) and columns (c).
An arb_mat_t is defined as an array of length one of type arb_mat_struct, permitting an arb_mat_t to be passed by reference.
-
arb_mat_entry(mat, i, j)¶
Macro giving a pointer to the entry at row i and column j.
-
arb_mat_nrows(mat)¶
Returns the number of rows of the matrix.
-
arb_mat_ncols(mat)¶
Returns the number of columns of the matrix.
Memory management¶
-
void arb_mat_init(arb_mat_t mat, slong r, slong c)¶
Initializes the matrix, setting it to the zero matrix with r rows and c columns.
-
slong arb_mat_allocated_bytes(const arb_mat_t x)¶
Returns the total number of bytes heap-allocated internally by this object. The count excludes the size of the structure itself. Add
sizeof(arb_mat_struct)to get the size of the object as a whole.
Conversions¶
-
void arb_mat_set_fmpz_mat(arb_mat_t dest, const fmpz_mat_t src)¶
-
void arb_mat_set_round_fmpz_mat(arb_mat_t dest, const fmpz_mat_t src, slong prec)¶
-
void arb_mat_set_fmpq_mat(arb_mat_t dest, const fmpq_mat_t src, slong prec)¶
Sets dest to src. The operands must have identical dimensions.
Random generation¶
Input and output¶
Comparisons¶
Predicate methods return 1 if the property certainly holds and 0 otherwise.
-
int arb_mat_equal(const arb_mat_t mat1, const arb_mat_t mat2)¶
Returns whether the matrices have the same dimensions and identical intervals as entries.
-
int arb_mat_overlaps(const arb_mat_t mat1, const arb_mat_t mat2)¶
Returns whether the matrices have the same dimensions and each entry in mat1 overlaps with the corresponding entry in mat2.
-
int arb_mat_contains_fmpz_mat(const arb_mat_t mat1, const fmpz_mat_t mat2)¶
-
int arb_mat_contains_fmpq_mat(const arb_mat_t mat1, const fmpq_mat_t mat2)¶
Returns whether the matrices have the same dimensions and each entry in mat2 is contained in the corresponding entry in mat1.
-
int arb_mat_eq(const arb_mat_t mat1, const arb_mat_t mat2)¶
Returns whether mat1 and mat2 certainly represent the same matrix.
-
int arb_mat_ne(const arb_mat_t mat1, const arb_mat_t mat2)¶
Returns whether mat1 and mat2 certainly do not represent the same matrix.
-
int arb_mat_is_empty(const arb_mat_t mat)¶
Returns whether the number of rows or the number of columns in mat is zero.
-
int arb_mat_is_square(const arb_mat_t mat)¶
Returns whether the number of rows is equal to the number of columns in mat.
-
int arb_mat_is_triu(const arb_mat_t mat)¶
Returns whether mat is upper triangular; that is, all entries below the main diagonal are exactly zero.
Special matrices¶
-
void arb_mat_one(arb_mat_t mat)¶
Sets the entries on the main diagonal to ones, and all other entries to zero.
-
void arb_mat_hilbert(arb_mat_t mat, slong prec)¶
Sets mat to the Hilbert matrix, which has entries \(A_{j,k} = 1/(j+k+1)\).
-
void arb_mat_pascal(arb_mat_t mat, int triangular, slong prec)¶
Sets mat to a Pascal matrix, whose entries are binomial coefficients. If triangular is 0, constructs a full symmetric matrix with the rows of Pascal’s triangle as successive antidiagonals. If triangular is 1, constructs the upper triangular matrix with the rows of Pascal’s triangle as columns, and if triangular is -1, constructs the lower triangular matrix with the rows of Pascal’s triangle as rows.
The entries are computed using recurrence relations. When the dimensions get large, some precision loss is possible; in that case, the user may wish to create the matrix at slightly higher precision and then round it to the final precision.
-
void arb_mat_stirling(arb_mat_t mat, int kind, slong prec)¶
Sets mat to a Stirling matrix, whose entries are Stirling numbers. If kind is 0, the entries are set to the unsigned Stirling numbers of the first kind. If kind is 1, the entries are set to the signed Stirling numbers of the first kind. If kind is 2, the entries are set to the Stirling numbers of the second kind.
The entries are computed using recurrence relations. When the dimensions get large, some precision loss is possible; in that case, the user may wish to create the matrix at slightly higher precision and then round it to the final precision.
-
void arb_mat_dct(arb_mat_t mat, int type, slong prec)¶
Sets mat to the DCT (discrete cosine transform) matrix of order n where n is the smallest dimension of mat (if mat is not square, the matrix is extended periodically along the larger dimension). There are many different conventions for defining DCT matrices; here, we use the normalized “DCT-II” transform matrix
\[A_{j,k} = \sqrt{\frac{2}{n}} \cos\left(\frac{\pi j}{n} \left(k+\frac{1}{2}\right)\right)\]which satisfies \(A^{-1} = A^T\). The type parameter is currently ignored and should be set to 0. In the future, it might be used to select a different convention.
Transpose¶
Norms¶
-
void arb_mat_bound_inf_norm(mag_t b, const arb_mat_t A)¶
Sets b to an upper bound for the infinity norm (i.e. the largest absolute value row sum) of A.
Arithmetic¶
-
void arb_mat_neg(arb_mat_t dest, const arb_mat_t src)¶
Sets dest to the exact negation of src. The operands must have the same dimensions.
-
void arb_mat_add(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec)¶
Sets res to the sum of mat1 and mat2. The operands must have the same dimensions.
-
void arb_mat_sub(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec)¶
Sets res to the difference of mat1 and mat2. The operands must have the same dimensions.
-
void arb_mat_mul(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec)¶
Sets res to the matrix product of mat1 and mat2. The operands must have compatible dimensions for matrix multiplication.
The classical version performs matrix multiplication in the trivial way.
The block version decomposes the input matrices into one or several blocks of uniformly scaled matrices and multiplies large blocks via fmpz_mat_mul. It also invokes
_arb_mat_addmul_rad_mag_fast()for the radius matrix multiplications.The threaded version performs classical multiplication but splits the computation over the number of threads returned by flint_get_num_threads().
The default version chooses an algorithm automatically.
-
void arb_mat_mul_entrywise(arb_mat_t C, const arb_mat_t A, const arb_mat_t B, slong prec)¶
Sets C to the entrywise product of A and B. The operands must have the same dimensions.
-
void arb_mat_sqr(arb_mat_t res, const arb_mat_t mat, slong prec)¶
Sets res to the matrix square of mat. The operands must both be square with the same dimensions.
-
void arb_mat_pow_ui(arb_mat_t res, const arb_mat_t mat, ulong exp, slong prec)¶
Sets res to mat raised to the power exp. Requires that mat is a square matrix.
-
void _arb_mat_addmul_rad_mag_fast(arb_mat_t C, mag_srcptr A, mag_srcptr B, slong ar, slong ac, slong bc)¶
Helper function for matrix multiplication. Adds to the radii of C the matrix product of the matrices represented by A and B, where A is a linear array of coefficients in row-major order and B is a linear array of coefficients in column-major order. This function assumes that all exponents are small and is unsafe for general use.
Scalar arithmetic¶
-
void arb_mat_scalar_mul_2exp_si(arb_mat_t B, const arb_mat_t A, slong c)¶
Sets B to A multiplied by \(2^c\).
-
void arb_mat_scalar_addmul_arb(arb_mat_t B, const arb_mat_t A, const arb_t c, slong prec)¶
Sets B to \(B + A \times c\).
Gaussian elimination and solving¶
-
int arb_mat_lu(slong *perm, arb_mat_t LU, const arb_mat_t A, slong prec)¶
Given an \(n \times n\) matrix \(A\), computes an LU decomposition \(PLU = A\) using Gaussian elimination with partial pivoting. The input and output matrices can be the same, performing the decomposition in-place.
Entry \(i\) in the permutation vector perm is set to the row index in the input matrix corresponding to row \(i\) in the output matrix.
The algorithm succeeds and returns nonzero if it can find \(n\) invertible (i.e. not containing zero) pivot entries. This guarantees that the matrix is invertible.
The algorithm fails and returns zero, leaving the entries in \(P\) and \(LU\) undefined, if it cannot find \(n\) invertible pivot elements. In this case, either the matrix is singular, the input matrix was computed to insufficient precision, or the LU decomposition was attempted at insufficient precision.
The classical version uses Gaussian elimination directly while the recursive version performs the computation in a block recursive way to benefit from fast matrix multiplication. The default version chooses an algorithm automatically.
-
void arb_mat_solve_tril_classical(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, int unit, slong prec)¶
-
void arb_mat_solve_tril_recursive(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, int unit, slong prec)¶
-
void arb_mat_solve_triu_classical(arb_mat_t X, const arb_mat_t U, const arb_mat_t B, int unit, slong prec)¶
-
void arb_mat_solve_triu_recursive(arb_mat_t X, const arb_mat_t U, const arb_mat_t B, int unit, slong prec)¶
-
void arb_mat_solve_triu(arb_mat_t X, const arb_mat_t U, const arb_mat_t B, int unit, slong prec)¶
Solves the lower triangular system \(LX = B\) or the upper triangular system \(UX = B\), respectively. If unit is set, the main diagonal of L or U is taken to consist of all ones, and in that case the actual entries on the diagonal are not read at all and can contain other data.
The classical versions perform the computations iteratively while the recursive versions perform the computations in a block recursive way to benefit from fast matrix multiplication. The default versions choose an algorithm automatically.
-
void arb_mat_solve_lu_precomp(arb_mat_t X, const slong *perm, const arb_mat_t LU, const arb_mat_t B, slong prec)¶
Solves \(AX = B\) given the precomputed nonsingular LU decomposition \(A = PLU\). The matrices \(X\) and \(B\) are allowed to be aliased with each other, but \(X\) is not allowed to be aliased with \(LU\).
-
int arb_mat_solve_precond(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec)¶
Solves \(AX = B\) where \(A\) is a nonsingular \(n \times n\) matrix and \(X\) and \(B\) are \(n \times m\) matrices.
If \(m > 0\) and \(A\) cannot be inverted numerically (indicating either that \(A\) is singular or that the precision is insufficient), the values in the output matrix are left undefined and zero is returned. A nonzero return value guarantees that \(A\) is invertible and that the exact solution matrix is contained in the output.
Three algorithms are provided:
The lu version performs LU decomposition directly in ball arithmetic. This is fast, but the bounds typically blow up exponentially with n, even if the system is well-conditioned. This algorithm is usually the best choice at very high precision.
The precond version computes an approximate inverse to precondition the system [HS1967]. This is usually several times slower than direct LU decomposition, but the bounds do not blow up with n if the system is well-conditioned. This algorithm is usually the best choice for large systems at low to moderate precision.
The default version selects between lu and precomp automatically.
The automatic choice should be reasonable most of the time, but users may benefit from trying either lu or precond in specific applications. For example, the lu solver often performs better for ill-conditioned systems where use of very high precision is unavoidable.
-
int arb_mat_solve_preapprox(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, const arb_mat_t R, const arb_mat_t T, slong prec)¶
Solves \(AX = B\) where \(A\) is a nonsingular \(n \times n\) matrix and \(X\) and \(B\) are \(n \times m\) matrices, given an approximation \(R\) of the matrix inverse of \(A\), and given the approximation \(T\) of the solution \(X\).
If \(m > 0\) and \(A\) cannot be inverted numerically (indicating either that \(A\) is singular or that the precision is insufficient, or that \(R\) is not a close enough approximation of the inverse of \(A\)), the values in the output matrix are left undefined and zero is returned. A nonzero return value guarantees that \(A\) is invertible and that the exact solution matrix is contained in the output.
-
int arb_mat_inv(arb_mat_t X, const arb_mat_t A, slong prec)¶
Sets \(X = A^{-1}\) where \(A\) is a square matrix, computed by solving the system \(AX = I\).
If \(A\) cannot be inverted numerically (indicating either that \(A\) is singular or that the precision is insufficient), the values in the output matrix are left undefined and zero is returned. A nonzero return value guarantees that the matrix is invertible and that the exact inverse is contained in the output.
-
void arb_mat_det(arb_t det, const arb_mat_t A, slong prec)¶
Sets det to the determinant of the matrix A.
The lu version uses Gaussian elimination with partial pivoting. If at some point an invertible pivot element cannot be found, the elimination is stopped and the magnitude of the determinant of the remaining submatrix is bounded using Hadamard’s inequality.
The precond version computes an approximate LU factorization of A and multiplies by the inverse L and U martices as preconditioners to obtain a matrix close to the identity matrix [Rum2010]. An enclosure for this determinant is computed using Gershgorin circles. This is about four times slower than direct Gaussian elimination, but much more numerically stable.
The default version automatically selects between the lu and precond versions and additionally handles small or triangular matrices by direct formulas.
-
void arb_mat_approx_solve_triu(arb_mat_t X, const arb_mat_t U, const arb_mat_t B, int unit, slong prec)¶
-
void arb_mat_approx_solve_tril(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, int unit, slong prec)¶
-
void arb_mat_approx_solve_lu_precomp(arb_mat_t X, const slong *perm, const arb_mat_t A, const arb_mat_t B, slong prec)¶
-
int arb_mat_approx_inv(arb_mat_t X, const arb_mat_t A, slong prec)¶
These methods perform approximate solving without any error control. The radii in the input matrices are ignored, the computations are done numerically with floating-point arithmetic (using ordinary Gaussian elimination and triangular solving, accelerated through the use of block recursive strategies for large matrices), and the output matrices are set to the approximate floating-point results with zeroed error bounds.
Approximate solutions are useful for computing preconditioning matrices for certified solutions. Some users may also find these methods useful for doing ordinary numerical linear algebra in applications where error bounds are not needed.
Cholesky decomposition and solving¶
-
int arb_mat_cho(arb_mat_t L, const arb_mat_t A, slong prec)¶
Computes the Cholesky decomposition of A, returning nonzero iff the symmetric matrix defined by the lower triangular part of A is certainly positive definite.
If a nonzero value is returned, then L is set to the lower triangular matrix such that \(A = L * L^T\).
If zero is returned, then either the matrix is not symmetric positive definite, the input matrix was computed to insufficient precision, or the decomposition was attempted at insufficient precision.
The underscore method computes L from A in-place, leaving the strict upper triangular region undefined.
-
void arb_mat_solve_cho_precomp(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, slong prec)¶
Solves \(AX = B\) given the precomputed Cholesky decomposition \(A = L L^T\). The matrices X and B are allowed to be aliased with each other, but X is not allowed to be aliased with L.
-
int arb_mat_spd_solve(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec)¶
Solves \(AX = B\) where A is a symmetric positive definite matrix and X and B are \(n \times m\) matrices, using Cholesky decomposition.
If \(m > 0\) and A cannot be factored using Cholesky decomposition (indicating either that A is not symmetric positive definite or that the precision is insufficient), the values in the output matrix are left undefined and zero is returned. A nonzero return value guarantees that the symmetric matrix defined through the lower triangular part of A is invertible and that the exact solution matrix is contained in the output.
-
void arb_mat_inv_cho_precomp(arb_mat_t X, const arb_mat_t L, slong prec)¶
Sets \(X = A^{-1}\) where \(A\) is a symmetric positive definite matrix whose Cholesky decomposition L has been computed with
arb_mat_cho(). The inverse is calculated using the method of [Kri2013] which is more efficient than solving \(AX = I\) witharb_mat_solve_cho_precomp().
-
int arb_mat_spd_inv(arb_mat_t X, const arb_mat_t A, slong prec)¶
Sets \(X = A^{-1}\) where A is a symmetric positive definite matrix. It is calculated using the method of [Kri2013] which computes fewer intermediate results than solving \(AX = I\) with
arb_mat_spd_solve().If A cannot be factored using Cholesky decomposition (indicating either that A is not symmetric positive definite or that the precision is insufficient), the values in the output matrix are left undefined and zero is returned. A nonzero return value guarantees that the symmetric matrix defined through the lower triangular part of A is invertible and that the exact inverse is contained in the output.
-
int arb_mat_ldl(arb_mat_t res, const arb_mat_t A, slong prec)¶
Computes the \(LDL^T\) decomposition of A, returning nonzero iff the symmetric matrix defined by the lower triangular part of A is certainly positive definite.
If a nonzero value is returned, then res is set to a lower triangular matrix that encodes the \(L * D * L^T\) decomposition of A. In particular, \(L\) is a lower triangular matrix with ones on its diagonal and whose strictly lower triangular region is the same as that of res. \(D\) is a diagonal matrix with the same diagonal as that of res.
If zero is returned, then either the matrix is not symmetric positive definite, the input matrix was computed to insufficient precision, or the decomposition was attempted at insufficient precision.
The underscore methods compute res from A in-place, leaving the strict upper triangular region undefined. The default method uses algorithm 4.1.2 from [GVL1996].
-
void arb_mat_solve_ldl_precomp(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, slong prec)¶
Solves \(AX = B\) given the precomputed \(A = LDL^T\) decomposition encoded by L. The matrices X and B are allowed to be aliased with each other, but X is not allowed to be aliased with L.
-
void arb_mat_inv_ldl_precomp(arb_mat_t X, const arb_mat_t L, slong prec)¶
Sets \(X = A^{-1}\) where \(A\) is a symmetric positive definite matrix whose \(LDL^T\) decomposition encoded by L has been computed with
arb_mat_ldl(). The inverse is calculated using the method of [Kri2013] which is more efficient than solving \(AX = I\) witharb_mat_solve_ldl_precomp().
Characteristic polynomial and companion matrix¶
-
void arb_mat_charpoly(arb_poly_t poly, const arb_mat_t mat, slong prec)¶
Sets poly to the characteristic polynomial of mat which must be a square matrix. If the matrix has n rows, the underscore method requires space for \(n + 1\) output coefficients. Employs a division-free algorithm using \(O(n^4)\) operations.
-
void _arb_mat_companion(arb_mat_t mat, arb_srcptr poly, slong prec)¶
-
void arb_mat_companion(arb_mat_t mat, const arb_poly_t poly, slong prec)¶
Sets the n by n matrix mat to the companion matrix of the polynomial poly which must have degree n. The underscore method reads \(n + 1\) input coefficients.
Special functions¶
-
void arb_mat_exp_taylor_sum(arb_mat_t S, const arb_mat_t A, slong N, slong prec)¶
Sets S to the truncated exponential Taylor series \(S = \sum_{k=0}^{N-1} A^k / k!\). Uses rectangular splitting to compute the sum using \(O(\sqrt{N})\) matrix multiplications. The recurrence relation for factorials is used to get scalars that are small integers instead of full factorials. As in [Joh2014b], all divisions are postponed to the end by computing partial factorials of length \(O(\sqrt{N})\). The scalars could be reduced by doing more divisions, but this appears to be slower in most cases.
-
void arb_mat_exp(arb_mat_t B, const arb_mat_t A, slong prec)¶
Sets B to the exponential of the matrix A, defined by the Taylor series
\[\exp(A) = \sum_{k=0}^{\infty} \frac{A^k}{k!}.\]The function is evaluated as \(\exp(A/2^r)^{2^r}\), where \(r\) is chosen to give rapid convergence.
The elementwise error when truncating the Taylor series after N terms is bounded by the error in the infinity norm, for which we have
\[\left\|\exp(2^{-r}A) - \sum_{k=0}^{N-1} \frac{\left(2^{-r} A\right)^k}{k!} \right\|_{\infty} = \left\|\sum_{k=N}^{\infty} \frac{\left(2^{-r} A\right)^k}{k!}\right\|_{\infty} \le \sum_{k=N}^{\infty} \frac{(2^{-r} \|A\|_{\infty})^k}{k!}.\]We bound the sum on the right using
mag_exp_tail(). Truncation error is not added to entries whose values are determined by the sparsity structure of \(A\).
Sparsity structure¶
-
void arb_mat_entrywise_is_zero(fmpz_mat_t dest, const arb_mat_t src)¶
Sets each entry of dest to indicate whether the corresponding entry of src is certainly zero. If the entry of src at row \(i\) and column \(j\) is zero according to
arb_is_zero()then the entry of dest at that row and column is set to one, otherwise that entry of dest is set to zero.
-
void arb_mat_entrywise_not_is_zero(fmpz_mat_t dest, const arb_mat_t src)¶
Sets each entry of dest to indicate whether the corresponding entry of src is not certainly zero. This the complement of
arb_mat_entrywise_is_zero().
-
slong arb_mat_count_is_zero(const arb_mat_t mat)¶
Returns the number of entries of mat that are certainly zero according to
arb_is_zero().
Component and error operations¶
Eigenvalues and eigenvectors¶
To compute eigenvalues and eigenvectors, one can convert to an
acb_mat_t and use the functions in acb_mat.h: Eigenvalues and eigenvectors.
In the future dedicated methods for real matrices will be added here.