Linear Algebra
np_matmul (a, b) — Function
Matrix multiplication. Supports complex arrays.
Computes the matrix product of two 2D ndarrays using BLAS. The number of columns of a must equal the number of rows of b.
Examples
(%i1) A : ndarray(matrix([1, 2], [3, 4]));
(%o1) ndarray([2, 2], DOUBLE-FLOAT)
(%i2) I : np_eye(2);
(%o2) ndarray([2, 2], DOUBLE-FLOAT)
(%i3) np_to_matrix(np_matmul(I, A));
(%o3) matrix([1.0, 2.0], [3.0, 4.0])
(%i4) np_to_matrix(np_matmul(A, A));
(%o4) matrix([7.0, 10.0], [15.0, 22.0])
See also: np_dot, np_inv, np_solve
np_inv (a) — Function
Matrix inverse. Supports complex arrays.
Computes the inverse of a square matrix using LAPACK. Signals an error if the matrix is singular or nearly singular.
Examples
(%i1) A : ndarray(matrix([1, 2], [3, 4]));
(%o1) ndarray([2, 2], DOUBLE-FLOAT)
(%i2) Ainv : np_inv(A);
(%o2) ndarray([2, 2], DOUBLE-FLOAT)
(%i3) np_to_matrix(np_matmul(A, Ainv));
(%o3) matrix([1.0, 0.0], [0.0, 1.0])
(%i4) errcatch(np_inv(ndarray(matrix([1, 0], [0, 0]))));
(%o4) []
See also: np_solve, np_det, np_pinv
np_det (a) — Function
Determinant of a square matrix.
Returns a scalar (complex for complex input). Uses LAPACK for computation.
Examples
(%i1) A : ndarray(matrix([1, 2], [3, 4]));
(%o1) ndarray([2, 2], DOUBLE-FLOAT)
(%i2) np_det(A);
(%o2) -2.0
(%i3) np_det(np_eye(5));
(%o3) 1.0
See also: np_inv, np_rank
np_solve (a, b) — Function
Solve the linear system Ax = b. Supports complex arrays.
Computes the solution to a system of linear equations where a is a square coefficient matrix and b is a right-hand side vector or matrix. Uses LAPACK.
Examples
(%i1) A : ndarray(matrix([2, 1], [5, 3]));
(%o1) ndarray([2, 2], DOUBLE-FLOAT)
(%i2) b : ndarray(matrix([4], [7]));
(%o2) ndarray([2, 1], DOUBLE-FLOAT)
(%i3) x : np_solve(A, b);
(%o3) ndarray([2, 1], DOUBLE-FLOAT)
(%i4) np_to_matrix(np_matmul(A, x));
(%o4) matrix([4.0], [7.0])
See also: np_inv, np_lstsq
np_svd (a) — Function
Singular Value Decomposition. Supports complex arrays.
Decomposes a into U, S, and Vt such that A = U * diag(S) * Vt, where U and Vt are orthogonal (or unitary for complex) matrices. S is returned as a 1D ndarray of singular values (not a diagonal matrix), always double-float. Returns a Maxima list [U, S, Vt].
Returns the economy (reduced) SVD: for an m-by-n matrix with k = min(m,n), U is m-by-k, S has k elements, and Vt is k-by-n. This means np_matmul(np_matmul(U, np_diag(np_to_list(S))), Vt) reconstructs A directly for any shape.
Examples
(%i1) A : ndarray(matrix([1, 0], [0, 2], [0, 0]));
(%o1) ndarray([3, 2], DOUBLE-FLOAT)
(%i2) [U, S, Vt] : np_svd(A);
(%o2) [ndarray, ndarray, ndarray]
(%i3) np_shape(U);
(%o3) [3, 2]
(%i4) np_shape(Vt);
(%o4) [2, 2]
(%i5) np_to_list(S);
(%o5) [2.0, 1.0]
(%i6) /* Reconstruct: U * diag(S) * Vt works directly */
np_norm(np_sub(A, np_matmul(np_matmul(U, np_diag(np_to_list(S))), Vt)));
(%o6) 0.0
See also: np_eig, np_rank, np_lstsq, np_pinv
np_eig (a) — Function
Eigendecomposition.
Computes eigenvalues and eigenvectors of a square matrix. Returns a Maxima list [eigenvalues, eigenvectors] where eigenvalues is a 1D ndarray and eigenvectors is a 2D ndarray whose columns are the eigenvectors.
When eigenvalues have non-negligible imaginary parts (e.g. rotation matrices, non-symmetric matrices), the returned ndarrays use complex-double-float dtype. When all eigenvalues are real, returns double-float ndarrays.
Examples
(%i1) A : ndarray(matrix([2, 1, 0], [1, 3, 1], [0, 1, 2]));
(%o1) ndarray([3, 3], DOUBLE-FLOAT)
(%i2) [vals, vecs] : np_eig(A);
(%o2) [ndarray, ndarray]
(%i3) np_to_list(vals);
(%o3) [1.0, 2.0, 4.0]
(%i4) /* Verify A*v = lambda*v */
np_to_matrix(np_sub(np_matmul(A, vecs),
np_matmul(vecs, np_diag(np_to_list(vals)))));
(%o4) matrix([0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0])
See also: np_svd
np_qr (a) — Function
QR decomposition. Works for any m-by-n matrix. Supports complex arrays.
Decomposes a into an orthogonal (or unitary for complex) matrix Q and an upper triangular (or upper trapezoidal) matrix R such that A = Q * R. Returns a Maxima list [Q, R]. For tall/square matrices (m >= n), uses LAPACK. For wide matrices (m < n), uses Householder reflections.
Examples
(%i1) A : ndarray(matrix([1, 2], [3, 4], [5, 6]));
(%o1) ndarray([3, 2], DOUBLE-FLOAT)
(%i2) [Q, R] : np_qr(A);
(%o2) [ndarray, ndarray]
(%i3) np_shape(Q);
(%o3) [3, 3]
(%i4) np_shape(R);
(%o4) [3, 2]
See also: np_lu, np_svd
np_lu (a) — Function
LU decomposition with partial pivoting. Supports complex arrays.
Decomposes a into a lower triangular matrix L (unit diagonal), an upper triangular matrix U, and a permutation matrix P such that P * A = L * U. Returns a Maxima list [L, U, P].
For an m-by-n matrix, L is m-by-min(m,n), U is min(m,n)-by-n, and P is m-by-m.
Examples
(%i1) A : ndarray(matrix([1, 2, 3], [4, 5, 6], [7, 8, 10]));
(%o1) ndarray([3, 3], DOUBLE-FLOAT)
(%i2) [L, U, P] : np_lu(A);
(%o2) [ndarray, ndarray, ndarray]
(%i3) /* Verify P*A = L*U */
np_to_matrix(np_sub(np_matmul(P, A), np_matmul(L, U)));
(%o3) matrix([0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0])
See also: np_qr, np_svd, np_solve
np_norm (a) / np_norm (a, ord) — Function
Matrix or vector norm. Supports complex arrays.
By default, computes the 2-norm for vectors and the Frobenius norm for matrices. The optional ord parameter selects the norm type. The result is always real (double-float).
Calling forms:
np_norm(a)– default norm (2-norm for vectors, Frobenius for matrices)np_norm(a, 1)– 1-norm (sum of abs for vectors, max column sum for matrices)np_norm(a, 2)– 2-norm (Euclidean for vectors, spectral/largest singular value for matrices)np_norm(a, inf)– infinity norm (max abs for vectors, max row sum for matrices)np_norm(a, fro)– Frobenius norm (matrices only, same as default)
Examples
(%i1) v : ndarray([1, -2, 3], [3]);
(%o1) ndarray([3], DOUBLE-FLOAT)
(%i2) np_norm(v);
(%o2) 3.7416573867739413
(%i3) np_norm(v, 1);
(%o3) 6.0
(%i4) np_norm(v, inf);
(%o4) 3.0
(%i5) A : ndarray(matrix([1, -7], [2, -3]));
(%o5) ndarray([2, 2], DOUBLE-FLOAT)
(%i6) np_norm(A, 1);
(%o6) 10.0
(%i7) np_norm(A, inf);
(%o7) 8.0
(%i8) np_norm(A, 2);
(%o8) 7.649700064568801
See also: np_det, np_rank
np_rank (a) — Function
Numerical rank of a matrix.
Computed via SVD. Counts the number of singular values above a tolerance threshold (proportional to machine epsilon and the largest singular value). Returns an integer.
Examples
(%i1) np_rank(np_eye(3));
(%o1) 3
(%i2) np_rank(ndarray(matrix([1, 2], [2, 4])));
(%o2) 1
(%i3) np_rank(np_zeros([3, 3]));
(%o3) 0
See also: np_svd, np_det
np_trace (a) — Function
Matrix trace (sum of diagonal elements).
Returns a scalar equal to the sum of the elements on the main diagonal (complex for complex input).
Examples
(%i1) A : ndarray(matrix([1, 2], [3, 4]));
(%o1) ndarray([2, 2], DOUBLE-FLOAT)
(%i2) np_trace(A);
(%o2) 5.0
(%i3) np_trace(np_eye(5));
(%o3) 5.0
See also: np_det, np_diag
np_transpose (a) — Function
Matrix transpose. Supports complex arrays (transpose only, no conjugation; see np_ctranspose).
Returns a new ndarray with rows and columns swapped.
Examples
(%i1) A : ndarray(matrix([1, 2], [3, 4]));
(%o1) ndarray([2, 2], DOUBLE-FLOAT)
(%i2) np_to_matrix(np_transpose(A));
(%o2) matrix([1.0, 3.0], [2.0, 4.0])
(%i3) np_shape(np_transpose(ndarray(matrix([1, 2, 3]))));
(%o3) [3, 1]
See also: np_conj, np_ctranspose, np_reshape
np_conj (a) — Function
Element-wise complex conjugate.
Returns a new ndarray where each element is the complex conjugate of the corresponding element in a. For real (double-float) arrays, this returns a copy of a unchanged. For complex arrays, conjugates each element (negates the imaginary part).
Examples
(%i1) A : ndarray(matrix([1+%i, 2-3*%i]), complex);
(%o1) ndarray([1, 2], COMPLEX-DOUBLE-FLOAT)
(%i2) np_to_list(np_conj(A));
(%o2) [1.0 - 1.0*%i, 2.0 + 3.0*%i]
(%i3) /* Real arrays are unchanged */
np_to_list(np_conj(ndarray([1, 2, 3], [3])));
(%o3) [1.0, 2.0, 3.0]
See also: np_ctranspose, np_transpose, np_real, np_imag
np_ctranspose (a) — Function
Conjugate transpose (Hermitian transpose).
For real matrices, this is the same as np_transpose. For complex matrices, it transposes and takes the complex conjugate of each element.
Examples
(%i1) A : ndarray(matrix([1, 2], [3, 4]));
(%o1) ndarray([2, 2], DOUBLE-FLOAT)
(%i2) np_to_matrix(np_ctranspose(A));
(%o2) matrix([1.0, 3.0], [2.0, 4.0])
See also: np_conj, np_transpose
np_real (a) — Function
Extract real parts element-wise.
Returns a new double-float ndarray containing the real part of each element. For real arrays, this is equivalent to np_copy.
Examples
(%i1) A : ndarray(matrix([1+2*%i, 3-4*%i]), complex);
(%o1) ndarray([1, 2], COMPLEX-DOUBLE-FLOAT)
(%i2) np_to_list(np_real(A));
(%o2) [1.0, 3.0]
See also: np_imag, np_angle, np_conj
np_imag (a) — Function
Extract imaginary parts element-wise.
Returns a new double-float ndarray containing the imaginary part of each element. For real arrays, returns all zeros.
Examples
(%i1) A : ndarray(matrix([1+2*%i, 3-4*%i]), complex);
(%o1) ndarray([1, 2], COMPLEX-DOUBLE-FLOAT)
(%i2) np_to_list(np_imag(A));
(%o2) [2.0, -4.0]
See also: np_real, np_angle, np_conj
np_angle (a) — Function
Element-wise phase angle.
Returns a new double-float ndarray containing the phase angle (argument) of each element in radians. For real positive numbers, returns 0; for real negative numbers, returns pi.
Examples
(%i1) A : ndarray(matrix([1, -1, %i]), complex);
(%o1) ndarray([1, 3], COMPLEX-DOUBLE-FLOAT)
(%i2) np_to_list(np_angle(A));
(%o2) [0.0, 3.141592653589793, 1.5707963267948966]
See also: np_real, np_imag, np_abs
np_expm (a) — Function
Matrix exponential. Supports complex arrays.
Computes the matrix exponential e^A, which is defined as the infinite series I + A + A^2/2! + A^3/3! + …. This is different from element-wise np_exp, which applies the scalar exponential to each element independently.
Uses the scaling-and-squaring method with adaptive diagonal Pade approximation (orders 3, 5, 7, 9, or 13 selected automatically based on the matrix 1-norm). The algorithm follows Higham (2005), “The Scaling and Squaring Method for the Matrix Exponential Revisited”.
Examples
(%i1) /* expm(0) = I */
np_to_matrix(np_expm(np_zeros([2, 2])));
(%o1) matrix([1.0, 0.0], [0.0, 1.0])
(%i2) /* expm(I) = e*I */
np_ref(np_expm(np_eye(2)), 0, 0);
(%o2) 2.718281828459045
(%i3) /* Rotation matrix from skew-symmetric */
R : ndarray(matrix([0, -1], [1, 0]));
(%o3) ndarray([2, 2], DOUBLE-FLOAT)
(%i4) E : np_to_matrix(np_expm(R));
(%o4) matrix([0.5403023058681398, -0.8414709848078965],
[0.8414709848078965, 0.5403023058681398])
(%i5) /* expm(A) * expm(-A) = I */
A : ndarray(matrix([1, 2], [3, 4]));
(%o5) ndarray([2, 2], DOUBLE-FLOAT)
(%i6) np_to_matrix(np_matmul(np_expm(A), np_expm(np_neg(A))));
(%o6) matrix([1.0, 0.0], [0.0, 1.0])
See also: np_exp, np_matmul
np_lstsq (a, b) — Function
Least-squares solution to Ax = b. Supports complex arrays.
Finds x that minimizes ||Ax - b||_2 using SVD. Works for over-determined systems (more equations than unknowns) and under-determined systems. For square full-rank matrices, gives the same result as np_solve.
Returns a Maxima list [x, residuals, rank, S] where:
- x — n-by-p solution ndarray
- residuals — 1D ndarray of squared residual norms
||Ax_j - b_j||^2for each column j of b (alwaysdouble-float). Only non-empty whenm > n(overdetermined) and A is full rank; otherwise an empty list[]. - rank — effective rank of A (integer)
- S — 1D ndarray of singular values of A (always
double-float)
Breaking change: Previous versions returned only x. Update callers from x : np_lstsq(A, b) to [x, residuals, rank, S] : np_lstsq(A, b).
Examples
(%i1) /* Linear fit: y = a + b*x through (1,1), (2,2), (3,3) */
A : ndarray(matrix([1, 1], [1, 2], [1, 3]));
(%o1) ndarray([3, 2], DOUBLE-FLOAT)
(%i2) b : ndarray(matrix([1], [2], [3]));
(%o2) ndarray([3, 1], DOUBLE-FLOAT)
(%i3) [x, residuals, rank, S] : np_lstsq(A, b);
(%o3) [ndarray, ndarray, 2, ndarray]
(%i4) np_to_list(x);
(%o4) [0.0, 1.0]
(%i5) rank;
(%o5) 2
(%i6) np_to_list(S);
(%o6) [3.86..., 0.64...]
See also: np_solve, np_pinv, np_svd
np_pinv (a) — Function
Moore-Penrose pseudo-inverse. Supports complex arrays.
Computed via SVD as A+ = V * S+ * Ut, where S+ inverts the non-zero singular values. For invertible square matrices, this is equivalent to np_inv. For non-square or rank-deficient matrices, it gives the best least-squares inverse.
The pseudo-inverse satisfies the four Moore-Penrose conditions: AA+A = A, A+AA+ = A+, (AA+)^T = AA+, and (A+*A)^T = A+*A.
For an m-by-n matrix, the pseudo-inverse is n-by-m.
Examples
(%i1) A : ndarray(matrix([1, 2], [3, 4], [5, 6]));
(%o1) ndarray([3, 2], DOUBLE-FLOAT)
(%i2) Ap : np_pinv(A);
(%o2) ndarray([2, 3], DOUBLE-FLOAT)
(%i3) /* Verify A * pinv(A) * A = A */
np_to_matrix(np_sub(np_matmul(np_matmul(A, Ap), A), A));
(%o3) matrix([0.0, 0.0], [0.0, 0.0], [0.0, 0.0])
(%i4) /* pinv(A) * b gives least-squares solution */
b : ndarray(matrix([1], [2], [3]));
(%o4) ndarray([3, 1], DOUBLE-FLOAT)
(%i5) np_to_list(np_matmul(Ap, b));
(%o5) [0.0, 1.0]
See also: np_inv, np_lstsq, np_svd
np_outer (a, b) — Function
Outer product of two 1D ndarrays.
Returns an m-by-n 2D ndarray where element (i,j) is a[i] * b[j]. Both arguments must be 1D ndarrays.
Examples
(%i1) np_to_matrix(np_outer(ndarray([1, 2, 3], [3]), ndarray([4, 5], [2])));
(%o1) matrix([4.0, 5.0], [8.0, 10.0], [12.0, 15.0])
(%i2) np_shape(np_outer(np_ones([4]), np_ones([3])));
(%o2) [4, 3]
(%i3) /* Outer product with self gives rank-1 matrix */
v : ndarray([1, 2], [2]);
(%o3) ndarray([2], DOUBLE-FLOAT)
(%i4) np_det(np_outer(v, v));
(%o4) 0.0
See also: np_matmul, np_dot