vtl.la

Linear algebra wrappers over VSL. Functions accept VTL tensors and return VTL tensors; many decompositions currently return f64 results even when the input has another numeric type. Check each function's signature when preserving dtype matters.

matmul supports complex128 vector, matrix, and broadcast-batch products via a pure-V kernel. Complex vector products follow NumPy's convention and do not conjugate either operand. General eig accepts real matrices and can return complex eigenpairs; complex-input decomposition and solve routines remain unsupported.

import vtl
import vtl.la

a := vtl.from_2d[f64]([[2.0, 0.0], [0.0, 4.0]])!
b := vtl.from_2d[f64]([[1.0, 2.0], [3.0, 4.0]])!
c := la.matmul[f64](a, b)!
contracted := la.tensordot[f64](a, b, 1)!
determinant := la.det[f64](a)!
main_diagonal := la.diag[f64](a, 0)!

Available operations include dot, matmul, tensordot, diag, batched det and inv, trace, matrix norm, vector vector_norm and its axis variants, outer, cross, solve, lstsq, qr, lu, cholesky, pinv, matrix_rank, matrix_rank_batch, slogdet, svd, svdvals, batched matrix_norm, matrix_power, cond, symmetric eigh/eigvalsh, and general eig/eigvals. Shape and tolerance requirements are checked by each function.

solve(a, b) accepts square matrices or stacks of square matrices. Leading batch dimensions follow NumPy broadcasting. A vector right-hand side has shape [n] (or [..., n] for vector batches); a matrix right-hand side has shape [..., n, nrhs]. The result keeps the broadcast batch shape and removes the last axis for vector right-hand sides. Singular systems return an error.

lstsq(a, b) accepts a vector or matrix right-hand side. It preserves vector output shape, reports the effective numerical rank, returns squared residuals for overdetermined full-column-rank systems, and an empty residual tensor for rank-deficient or underdetermined systems.

tensordot(a, b, axes) contracts the last axes dimensions of a with the first axes dimensions of b. Use tensordot_axes(a, b, axes_a, axes_b) to choose arbitrary axis pairs.

vector_norm_axes(t, ord, axes, keepdims) reduces a tuple of axes in one operation. It requires at least one axis; axes may be negative but must be unique. For example, vector_norm_axes(t, 2, [0, 2], false) reduces the first and last dimensions and keeps the middle dimension.

diag(input, offset) builds a matrix from a vector or extracts a diagonal vector from a matrix. Positive offsets select diagonals above the main one; negative offsets select diagonals below it. It accepts rank-1 and rank-2 inputs and returns a copy.

covariance_matrix(data, rowvar, ddof) and correlation_matrix(data, rowvar) accept rank-2 tensors. With rowvar: true, each row is a variable and columns are observations, matching NumPy's default. Set rowvar: false when variables are columns. Covariance uses ddof: 0 for population values and ddof: 1 for sample values. Non-positive degrees of freedom follow NumPy's divisor semantics, including NaN when the divisor is zero. Correlation returns NaN for constant variables.

More examples: linear algebra tutorial and advanced tutorial.

svd returns (u, s, vt) for real matrices and batches. Its default full_matrices: true matches NumPy's square-factor shapes; set it to false for reduced factors. Results use f64, and inputs must be finite. The default build uses a pure-V Jacobi algorithm. Add -d vsl_lapack_lapacke to use VSL's native LAPACKE backend for SVD when LAPACKE is installed.

fn cholesky[T](a &vtl.Tensor[T]) !&vtl.Tensor[f64]

cholesky computes Cholesky factorization of a symmetric positive-definite matrix. Returns lower-triangular L where A = L * L^T.

fn cond[T](input &vtl.Tensor[T], options CondOptions) !&vtl.Tensor[f64]

fn correlation_matrix[T](data &vtl.Tensor[T], rowvar bool) !&vtl.Tensor[f64]

correlation_matrix returns the Pearson correlation matrix. As with covariance_matrix, rowvar selects whether observations occupy rows or columns. Constant variables produce NaN correlations, matching NumPy.

fn covariance_matrix[T](data &vtl.Tensor[T], rowvar bool, ddof int) !&vtl.Tensor[f64]

covariance_matrix returns a variables-by-variables covariance matrix. Input is interpreted like NumPy's cov: rows are variables when rowvar is true, and columns are variables when rowvar is false. ddof is the delta degrees of freedom: use 0 for population covariance and 1 for sample.

fn cross[T](u &vtl.Tensor[T], v &vtl.Tensor[T]) !&vtl.Tensor[f64]

cross computes the 3D cross product of two vectors: u × v Both vectors must have length 3.

fn det[T](t &vtl.Tensor[T]) !&vtl.Tensor[f64]

det exposes this operation as part of the public API.

fn det_complex(input &vtl.Tensor[vcomplex.Complex]) !&vtl.Tensor[vcomplex.Complex]

det_complex computes determinants of complex128 square matrices. Leading batch dimensions are preserved; singular matrices have determinant 0.

fn diag[T](input &vtl.Tensor[T], offset int) !&vtl.Tensor[T]

diag constructs a diagonal matrix from a vector or extracts a diagonal from a matrix. A positive offset moves above the main diagonal; a negative offset moves below it, matching NumPy's diag convention. The result is a copy.

fn dot[T](a &vtl.Tensor[T], b &vtl.Tensor[T]) !&vtl.Tensor[f64]

dot exposes this operation as part of the public API.

fn eig[T](input &vtl.Tensor[T]) !(&vtl.Tensor[vcomplex.Complex], &vtl.Tensor[vcomplex.Complex])

eig computes the general eigenvalues and right eigenvectors of each trailing real square matrix. Complex conjugate pairs are returned as complex values; eigenvectors are stored in columns. The outputs match NumPy shapes: [..., N] and [..., N, N].

fn eigh[T](input &vtl.Tensor[T], options EighOptions) !(&vtl.Tensor[f64], &vtl.Tensor[f64])

eigh computes ascending eigenvalues and corresponding eigenvectors for each trailing real symmetric matrix. Eigenvectors are stored in columns, and the results have shapes [..., N] and [..., N, N]. Non-finite inputs and matrices that do not converge within 100 Jacobi sweeps return errors.

fn eigvals[T](input &vtl.Tensor[T]) !&vtl.Tensor[vcomplex.Complex]

eigvals returns the general eigenvalues of each trailing real square matrix without computing eigenvectors.

fn eigvalsh[T](input &vtl.Tensor[T], options EighOptions) !&vtl.Tensor[f64]

eigvalsh returns only the ascending eigenvalues of trailing real symmetric matrices. The UPLO selection matches EighOptions and skips eigenvector work.

fn inv[T](t &vtl.Tensor[T]) !&vtl.Tensor[f64]

inv exposes this operation as part of the public API.

fn inv_complex(input &vtl.Tensor[vcomplex.Complex]) !&vtl.Tensor[vcomplex.Complex]

inv_complex computes the inverse of complex128 square matrices. Leading batch dimensions are preserved; singular matrices return an error.

fn lstsq[T](a &vtl.Tensor[T], b &vtl.Tensor[T]) !(&vtl.Tensor[f64], &vtl.Tensor[f64], int, &vtl.Tensor[f64])

lstsq solves the linear least-squares problem min ||Ax - B||_2. Returns (x, residuals, rank, singular_values).

fn lstsq_complex(a &vtl.Tensor[vcomplex.Complex], b &vtl.Tensor[vcomplex.Complex]) !(&vtl.Tensor[vcomplex.Complex], &vtl.Tensor[f64], int, &vtl.Tensor[f64])

lstsq_complex solves the complex linear least-squares problem min ||Ax - B||_2. It accepts a 2D A and vector or matrix B, and returns the solution, squared residual sums, numerical rank, and singular values. The current implementation embeds the problem in real arithmetic, so its temporary matrices are twice as large along each matrix dimension.

fn lu[T](a &vtl.Tensor[T]) !(&vtl.Tensor[f64], &vtl.Tensor[f64], &vtl.Tensor[int])

lu computes LU decomposition with partial pivoting: PA = LU. Returns (L, U, piv) as 2D tensors. L shape: [m, min(m,n)], U shape: [min(m,n), n]

fn matmul[T](a &vtl.Tensor[T], b &vtl.Tensor[T]) !&vtl.Tensor[T]

matmul exposes this operation as part of the public API.

fn matrix_norm[T](input &vtl.Tensor[T], options MatrixNormOptions) !&vtl.Tensor[f64]

matrix_norm computes a matrix norm for each trailing matrix in an N-D tensor. Supported orders are 'fro', 'nuc', '1', '-1', '2', '-2', 'inf', and '-inf'. The default Frobenius result is returned as a one-element tensor for a single matrix; keepdims retains the final two unit dimensions.

fn matrix_power[T](input &vtl.Tensor[T], exponent int) !&vtl.Tensor[f64]

matrix_power raises each trailing square matrix to an integer power using exponentiation by squaring. Negative powers use a partial-pivot inverse. Results are f64 tensors with the same shape as the input.

fn matrix_rank[T](a &vtl.Tensor[T], tol f64) !int

matrix_rank returns the effective numerical rank of A.

fn matrix_rank_batch[T](input &vtl.Tensor[T], tol f64) !&vtl.Tensor[int]

matrix_rank_batch returns the numerical rank of each trailing matrix. A positive tol is an absolute threshold; non-positive tol selects NumPy's dtype-aware default threshold.

fn multi_dot[T](operands []&vtl.Tensor[T]) !&vtl.Tensor[T]

multi_dot multiplies a sequence of matrices using the parenthesization with the lowest estimated scalar multiplication count. Vectors are allowed only as the first or last operand, matching NumPy's linalg.multi_dot contract.

fn norm[T](t &vtl.Tensor[T], ord string) !&vtl.Tensor[f64]

norm returns the matrix norm of a tensor. ord: "F" (Frobenius, default), "1" (column sum), "I" (row sum / infinity).

fn outer[T](u &vtl.Tensor[T], v &vtl.Tensor[T]) !&vtl.Tensor[f64]

outer computes the outer product of two vectors. result[i,j] = u[i] * v[j] result shape: [u.len, v.len]

fn pinv[T](a &vtl.Tensor[T], tol f64) !&vtl.Tensor[f64]

pinv computes the Moore-Penrose pseudoinverse of A using SVD.

fn qr[T](a &vtl.Tensor[T]) !(&vtl.Tensor[f64], &vtl.Tensor[f64])

qr computes QR factorization of A. Returns (Q, R) where Q is orthonormal and R is upper triangular. Q shape: [m, min(m,n)], R shape: [min(m,n), n]

fn slogdet[T](input &vtl.Tensor[T]) !(&vtl.Tensor[f64], &vtl.Tensor[f64])

slogdet returns each matrix's determinant sign and the natural logarithm of its absolute determinant. The two results have the leading batch shape; for a single matrix they contain one value. Input matrices must be finite and square, and may be stacked across leading dimensions.

fn slogdet_complex(input &vtl.Tensor[vcomplex.Complex]) !(&vtl.Tensor[vcomplex.Complex], &vtl.Tensor[f64])

slogdet_complex returns the unit-complex determinant phase and the natural logarithm of its magnitude for each complex128 square matrix. Leading batch dimensions are preserved. Singular matrices return a zero phase and -Inf.

fn solve[T](a &vtl.Tensor[T], b &vtl.Tensor[T]) !&vtl.Tensor[f64]

solve solves A * X = B for stacks of square A matrices. Leading dimensions broadcast like NumPy; B may be a vector or a matrix of right-hand sides.

fn solve_complex(a &vtl.Tensor[vcomplex.Complex], b &vtl.Tensor[vcomplex.Complex]) !&vtl.Tensor[vcomplex.Complex]

solve_complex solves A * X = B for complex128 square matrix stacks. Leading batch dimensions broadcast like NumPy; B may be a vector or matrix. Partial pivoting is used for each system. Singular systems return an error.

fn svd[T](input &vtl.Tensor[T], options SVDOptions) !(&vtl.Tensor[f64], &vtl.Tensor[f64], &vtl.Tensor[f64])

svd computes the singular value decomposition of every trailing real matrix, returning U, descending singular values, and V transpose. Results use f64. Reduced factors have shapes [..., M, K], [..., K], [..., K, N], where K=min(M,N); full factors use [..., M, M] and [..., N, N].

fn svdvals[T](input &vtl.Tensor[T]) !&vtl.Tensor[f64]

svdvals returns the singular values of every trailing real matrix, descending within each matrix. The output shape is [..., min(M, N)]. Non-finite inputs and inputs below rank two return errors.

fn tensordot[T](a &vtl.Tensor[T], b &vtl.Tensor[T], axes int) !&vtl.Tensor[T]

tensordot contracts the last axes dimensions of a with the first axes dimensions of b. The contracted dimensions must match pairwise. The result contains the remaining dimensions of a followed by the remaining dimensions of b, preserving the input element type.

fn tensordot_axes[T](a &vtl.Tensor[T], b &vtl.Tensor[T], axes_a []int, axes_b []int) !&vtl.Tensor[T]

tensordot_axes contracts the specified axes of a and b pairwise. Negative axis numbers count from the end, as in NumPy. Uncontracted dimensions are ordered as the remaining axes of a followed by the remaining dimensions of b.

fn trace[T](t &vtl.Tensor[T]) !&vtl.Tensor[f64]

trace returns the sum of diagonal elements of a square matrix.

fn trace_axes[T](input &vtl.Tensor[T], options TraceAxesOptions) !&vtl.Tensor[f64]

trace_axes sums a diagonal of each matrix selected by axis1 and axis2. Negative axes count from the end. The default axes are 0 and 1, matching numpy.trace; an N-D input returns the remaining dimensions as a batch.

fn vector_norm[T](t &vtl.Tensor[T], ord f64) !&vtl.Tensor[f64]

vector_norm computes a p-norm over all tensor elements and returns a one element f64 tensor. Positive and negative finite p, zero, and +/-infinity are supported. Scaling keeps finite results stable for very large or small data.

fn vector_norm_axes[T](t &vtl.Tensor[T], ord f64, axes []int, keepdims bool) !&vtl.Tensor[f64]

vector_norm_axes computes a p-norm over one or more axes. Negative axes count from the end; axes must be unique. At least one axis is required. With keepdims, reduced dimensions remain with length one.

fn vector_norm_axis[T](t &vtl.Tensor[T], ord f64, axis int) !&vtl.Tensor[f64]

vector_norm_axis computes a p-norm independently along axis and removes the reduced dimension. Negative axes count from the end.

fn vector_norm_axis_keepdims[T](t &vtl.Tensor[T], ord f64, axis int) !&vtl.Tensor[f64]

vector_norm_axis_keepdims retains the reduced dimension with length one.

struct CondOptions {
pub:
	ord string = '2'
}

cond computes a matrix condition number for each trailing square matrix. Supported orders match NumPy's real matrix orders: 2, -2, 1, -1, inf, -inf, and fro. A single matrix returns a one-element tensor.

struct EighOptions {
pub:
	uplo string = 'L'
}

EighOptions selects the triangle used to construct each symmetric matrix. The default lower triangle matches NumPy's default UPLO.

struct MatrixNormOptions {
pub:
	ord      string = 'fro'
	keepdims bool
}

MatrixNormOptions selects a NumPy-compatible matrix norm and output shape.

struct SVDOptions {
pub:
	full_matrices bool = true
}

SVDOptions selects whether svd returns full orthogonal matrices or the reduced factors. The default matches NumPy's full_matrices=true behavior.

struct TraceAxesOptions {
pub:
	offset int
	axis1  int
	axis2  int = 1
}

TraceAxesOptions selects the diagonal and the two axes traced over.