pub struct Matrix<T> { /* private fields */ }Expand description
A matrix of rows × cols elements of the type T, held column after column in one vector.
§Examples
use a::math::matrix::Matrix;
let a = Matrix::from_rows(&[[1.0, 2.0], [3.0, 4.0]]);
let b = Matrix::identity(2);
assert_eq!(&a * &b, a);
assert_eq!(a[(1, 0)], 3.0);
assert_eq!(a.transpose(), Matrix::from_rows(&[[1.0, 3.0], [2.0, 4.0]]));Implementations§
Source§impl<T: Element> Matrix<T>
impl<T: Element> Matrix<T>
Sourcepub fn conjugate_transpose(&self) -> Self
pub fn conjugate_transpose(&self) -> Self
The conjugate transpose A*, the transpose with each element conjugated: the transpose
itself for real elements.
Sourcepub fn hadamard(&self, other: &Matrix<T>) -> Matrix<T>
pub fn hadamard(&self, other: &Matrix<T>) -> Matrix<T>
The product elementwise, Hadamard’s A ∘ B.
§Panics
Panics if the matrices are not of the same shape.
Sourcepub fn add_scaled_row(&mut self, target: usize, source: usize, k: &T)
pub fn add_scaled_row(&mut self, target: usize, source: usize, k: &T)
Sourcepub fn swap_columns(&mut self, a: usize, b: usize)
pub fn swap_columns(&mut self, a: usize, b: usize)
Sourcepub fn scale_column(&mut self, j: usize, k: &T)
pub fn scale_column(&mut self, j: usize, k: &T)
Sourcepub fn add_scaled_column(&mut self, target: usize, source: usize, k: &T)
pub fn add_scaled_column(&mut self, target: usize, source: usize, k: &T)
Source§impl<T: Scalar> Matrix<T>
impl<T: Scalar> Matrix<T>
Sourcepub fn cholesky(&self) -> Option<Cholesky<T>>
pub fn cholesky(&self) -> Option<Cholesky<T>>
The Cholesky decomposition A = L·L*, or None when the matrix is not positive
definite: a pivot is not positive, or not a number. Only the lower triangle is read, the
upper one taken to be its conjugate. About n³/3 operations, nearly all of them in
products of blocks.
§Panics
Panics if the matrix is not square.
Source§impl<T: Scalar> Matrix<T>
impl<T: Scalar> Matrix<T>
Sourcepub fn hermitian_eigen(&self) -> HermitianEigen<T>
pub fn hermitian_eigen(&self) -> HermitianEigen<T>
The eigenvalues and orthonormal eigenvectors of a Hermitian (for real elements,
symmetric) matrix: see HermitianEigen. A blocked reduction to a real symmetric
tridiagonal matrix by Householder reflections, its eigensystem by Cuppen’s divide and
conquer, then the reflections applied to it in blocks; up to about a hundred rows (80 to
160 by the type), where that takes longer, the reduction one column at a time, its
reflections multiplied out into Q, and the implicit QL algorithm’s rotations applied
to Q. Only the lower triangle is read, the upper one taken to be its conjugate. A
matrix with a NaN or an infinity gives NaNs.
§Panics
Panics if the matrix is not square.
Sourcepub fn hermitian_eigenvalues(&self) -> Vec<T::Real>
pub fn hermitian_eigenvalues(&self) -> Vec<T::Real>
The eigenvalues of a Hermitian (for real elements, symmetric) matrix, real and
ascending, without the eigenvectors: the same reduction, then the QL and QR algorithm
without square roots, from about a hundred rows less than half the time of
hermitian_eigen. Only the lower triangle is read.
§Panics
Panics if the matrix is not square.
Source§impl Matrix<Number>
impl Matrix<Number>
Sourcepub fn det(&self) -> Number
pub fn det(&self) -> Number
The determinant, exactly, by Bareiss’s fraction-free elimination: each step’s elements
are minors of the matrix, (a_ij·a_kk - a_ik·a_kj) / p with the previous pivot p
dividing exactly, so the numbers grow no longer than the determinant. About n³
products. The determinant of the 0 × 0 matrix is 1.
§Panics
Panics if the matrix is not square.
§Examples
use a::math::matrix::Matrix;
use a::math::numbers::Number;
let a = Matrix::from_rows(&[[2, -1, 0], [-1, 2, -1], [0, -1, 2]]).map(|&x| Number::from(x));
assert_eq!(a.det(), 4);Sourcepub fn rank(&self) -> usize
pub fn rank(&self) -> usize
The rank, exactly: the number of pivots of Bareiss’s fraction-free elimination.
Sourcepub fn characteristic_polynomial(&self) -> Vec<Number>
pub fn characteristic_polynomial(&self) -> Vec<Number>
The coefficients of the characteristic polynomial det(x·I - A), from xⁿ down,
exactly, by Berkowitz’s method: the polynomial of each leading block from the one before
by a Toeplitz matrix of 1, -a_kk and -R·A_k^j·C, without a division. About n⁴/2
products.
§Panics
Panics if the matrix is not square.
§Examples
use a::math::matrix::Matrix;
use a::math::numbers::Number;
let a = Matrix::from_rows(&[[2, 1], [1, 2]]).map(|&x| Number::from(x));
assert_eq!(a.characteristic_polynomial(), [1, -4, 3].map(Number::from));Sourcepub fn adjugate(&self) -> Matrix<Number>
pub fn adjugate(&self) -> Matrix<Number>
The adjugate adj(A), the transpose of the matrix of cofactors, with
A·adj(A) = adj(A)·A = det(A)·I, for singular matrices too: by Cayley and Hamilton,
(-1)ⁿ⁻¹·(Aⁿ⁻¹ + c₁·Aⁿ⁻² + … + cₙ₋₁·I) for the characteristic polynomial’s c.
§Panics
Panics if the matrix is not square.
Sourcepub fn inverse(&self) -> Option<(Matrix<Number>, Number)>
pub fn inverse(&self) -> Option<(Matrix<Number>, Number)>
The inverse as the adjugate and the determinant, A⁻¹ = adj(A) / det(A), exactly, or
None when the determinant is zero.
§Panics
Panics if the matrix is not square.
§Examples
use a::math::matrix::Matrix;
use a::math::numbers::Number;
let a = Matrix::from_rows(&[[2, 1], [5, 3]]).map(|&x| Number::from(x));
let (adj, det) = a.inverse().unwrap();
assert_eq!(det, 1);
assert_eq!(adj, Matrix::from_rows(&[[3, -1], [-5, 2]]).map(|&x| Number::from(x)));Sourcepub fn null_space(&self) -> Matrix<Number>
pub fn null_space(&self) -> Matrix<Number>
A basis of the null space {x : A·x = 0} over the rationals, as the columns of an
n × (n - rank) matrix of integers: one vector for each column without a pivot in the
fraction-free echelon form, its free coordinate positive and its elements without a
common factor.
§Examples
use a::math::matrix::Matrix;
use a::math::numbers::Number;
let a = Matrix::from_rows(&[[1, 2, 3], [2, 4, 6]]).map(|&x| Number::from(x));
let null = a.null_space();
assert_eq!(null.shape(), (3, 2));
assert!((&a * &null).as_slice().iter().all(|x| x.is_zero()));Source§impl<T: Scalar> Matrix<T>
impl<T: Scalar> Matrix<T>
Sourcepub fn exp(&self) -> Self
pub fn exp(&self) -> Self
The exponential e^A = I + A + A²/2! + …, by scaling and squaring: A / 2^s small
enough for a Padé approximant of degree 3 to 13, chosen by its 1-norm (7 in single
precision), then squared s times, as Higham’s expm (2005).
§Panics
Panics if the matrix is not square.
§Examples
use a::math::matrix::Matrix;
// The generator of rotations: e^A rotates by θ.
let t = 0.5f64;
let a = Matrix::from_rows(&[[0.0, -t], [t, 0.0]]);
let r = Matrix::from_rows(&[[t.cos(), -t.sin()], [t.sin(), t.cos()]]);
assert!((&a.exp() - &r).as_slice().iter().all(|x| x.abs() < 1e-15));Sourcepub fn sqrt(&self) -> Option<Self>
pub fn sqrt(&self) -> Option<Self>
The principal square root √A, whose eigenvalues have positive real parts, by the
Schur method: R² = T solved for an upper triangular R element by element, then
Q·R·Q*. None for a real matrix with a negative real eigenvalue, whose square roots
are not real, or for a singular matrix that has none.
§Panics
Panics if the matrix is not square.
§Examples
use a::math::matrix::Matrix;
let a: Matrix<f64> = Matrix::from_rows(&[[33.0, 24.0], [48.0, 57.0]]);
let r = a.sqrt().unwrap();
assert!((&r - &Matrix::from_rows(&[[5.0, 2.0], [4.0, 7.0]])).as_slice().iter().all(|x| x.abs() < 1e-13));
assert!(Matrix::from_rows(&[[-1.0, 0.0], [0.0, 1.0]]).sqrt().is_none());Sourcepub fn log(&self) -> Option<Self>
pub fn log(&self) -> Option<Self>
The principal logarithm log A, whose eigenvalues have imaginary parts in (-π, π], by
inverse scaling and squaring: square roots of the Schur form T until it is near I,
then a Padé approximant of degree 7 in partial fractions, times 2^k for k square
roots, as Higham’s logm (2001) with Al-Mohy and Higham’s accurate diagonal. None
for a singular matrix, and for a real one with an eigenvalue on the negative real axis,
whose logarithms are not real.
§Panics
Panics if the matrix is not square.
§Examples
use a::math::matrix::Matrix;
let a: Matrix<f64> = Matrix::from_rows(&[[1.0, 2.0], [0.5, 3.0]]);
let back = a.log().unwrap().exp();
assert!((&back - &a).as_slice().iter().all(|x| x.abs() < 1e-13));
assert!(Matrix::<f64>::zeros(2, 2).log().is_none());Source§impl<T: Scalar> Matrix<T>
impl<T: Scalar> Matrix<T>
Sourcepub fn lu(&self) -> Lu<T>
pub fn lu(&self) -> Lu<T>
The LU decomposition with partial pivoting, P·A = L·U: see Lu. It takes about
2n³/3 operations for an n × n matrix, nearly all of them in products of blocks.
Sourcepub fn det(&self) -> T
pub fn det(&self) -> T
The determinant det(A), from the LU decomposition: the product of the pivots, with the
sign of the permutation. The determinant of the 0 × 0 matrix is 1.
§Panics
Panics if the matrix is not square.
Sourcepub fn inverse(&self) -> Option<Self>
pub fn inverse(&self) -> Option<Self>
The inverse A⁻¹, from the LU decomposition, or None when a pivot is zero: the matrix
is singular. A matrix that is singular only up to rounding gives an inverse with huge
elements; its condition number tells how far to trust it.
§Panics
Panics if the matrix is not square.
Source§impl<T: Scalar> Matrix<T>
impl<T: Scalar> Matrix<T>
Sourcepub fn norm(&self, norm: Norm) -> T::Real
pub fn norm(&self, norm: Norm) -> T::Real
The norm ‖A‖ of the kind given: see Norm. The 2-norm and the nuclear norm take the
singular values; the others one pass over the elements, Frobenius’s scaled so that the
squares neither overflow nor underflow.
§Examples
use a::math::matrix::{Matrix, Norm};
let a: Matrix<f64> = Matrix::from_rows(&[[1.0, -2.0], [-3.0, 4.0]]);
assert_eq!(a.norm(Norm::One), 6.0);
assert_eq!(a.norm(Norm::Infinity), 7.0);
assert_eq!(a.norm(Norm::Max), 4.0);
assert!((a.norm(Norm::Frobenius) - 30f64.sqrt()).abs() < 1e-15);Sourcepub fn cond(&self, norm: Norm) -> T::Real
pub fn cond(&self, norm: Norm) -> T::Real
The condition number κ(A) = ‖A‖·‖A⁻¹‖ in the norm given, which bounds how much a
relative change of A or of b can change the solution of A·x = b: infinite for a
singular matrix. In the 2-norm it is the ratio of the largest singular value to the
smallest, for a matrix of any shape; in the other norms it takes the inverse.
§Panics
Panics if the norm is not Norm::Two and the matrix is not square.
Sourcepub fn rank(&self) -> usize
pub fn rank(&self) -> usize
The rank: the number of singular values above max(m, n)·ε·σ₁, ε the unit roundoff
of the type and σ₁ the largest singular value, as NumPy’s matrix_rank.
Sourcepub fn null_space(&self) -> Matrix<T>
pub fn null_space(&self) -> Matrix<T>
An orthonormal basis of the null space {x : A·x = 0}, as the columns of an
n × (n - rank) matrix: the right singular vectors of the singular values at or below
the tolerance of rank.
Sourcepub fn characteristic_polynomial(&self) -> Vec<T>
pub fn characteristic_polynomial(&self) -> Vec<T>
The coefficients of the characteristic polynomial det(x·I - A), from xⁿ down: 1,
then -tr(A), …, (-1)ⁿ·det(A). From the Hessenberg form, by La Budde’s recurrence on
its leading blocks, about n³ operations.
§Panics
Panics if the matrix is not square.
§Examples
use a::math::matrix::Matrix;
let a: Matrix<f64> = Matrix::from_rows(&[[2.0, 1.0], [1.0, 2.0]]);
// (x - 1)(x - 3)
assert_eq!(a.characteristic_polynomial(), [1.0, -4.0, 3.0]);Source§impl<T: Scalar> Matrix<T>
impl<T: Scalar> Matrix<T>
Sourcepub fn schur(&self) -> Schur<T>
pub fn schur(&self) -> Schur<T>
The Schur decomposition A = Q·T·Q*: see Schur. A reduction to Hessenberg form by
blocks of Householder reflections, then the QR algorithm with implicit shifts, for large
matrices many at a time with aggressive early deflation, about 25n³ operations. A
matrix with a NaN or an infinity gives NaNs.
§Panics
Panics if the matrix is not square.
Sourcepub fn eigenvalues(&self) -> Vec<Complex<T::Real>>
pub fn eigenvalues(&self) -> Vec<Complex<T::Real>>
The eigenvalues, in no particular order, those of a real matrix that are not real in
pairs of conjugates: the diagonal of the Schur form, computed without Q, so about
10n³ operations. A matrix with a NaN or an infinity gives NaNs.
§Panics
Panics if the matrix is not square.
Source§impl<T: Scalar> Matrix<T>
impl<T: Scalar> Matrix<T>
Sourcepub fn svd(&self) -> Svd<T>
pub fn svd(&self) -> Svd<T>
The singular value decomposition A = U·Σ·V*: see Svd. A reduction to a real
bidiagonal matrix by Householder reflections from both sides, in panels whose updates of
the rest of the matrix are products of blocks, then divide and conquer on the bidiagonal
matrix, Gu and Eisenstat’s, and the reflections applied to its singular vectors in blocks:
about 9n³ operations for a square matrix, most of them in products of blocks. Up to
about a hundred columns (88 to 112 by the type), where that takes longer, the reduction
is made one column and row at a time, U and V are multiplied out from the reflections,
and the implicit QR algorithm’s rotations are applied to them. A matrix with at least
2.5 times as many rows as columns is first reduced to triangular by QR. A matrix with a
NaN or an infinity gives NaNs.
Sourcepub fn singular_values(&self) -> Vec<T::Real>
pub fn singular_values(&self) -> Vec<T::Real>
The singular values, real, not negative, in descending order, without U and V:
much less work than svd. Those of the bidiagonal matrix come from the
dqds algorithm, to high relative accuracy, so that a bidiagonal matrix has even its
smallest singular values to nearly full precision; QR goes first from 1.75 times as
many rows as columns.
Source§impl<T> Matrix<T>
impl<T> Matrix<T>
Sourcepub fn from_column_major(rows: usize, cols: usize, data: Vec<T>) -> Self
pub fn from_column_major(rows: usize, cols: usize, data: Vec<T>) -> Self
The matrix of rows × cols elements given column after column.
§Panics
Panics if data doesn’t hold rows · cols elements.
Sourcepub fn from_row_major(rows: usize, cols: usize, data: Vec<T>) -> Self
pub fn from_row_major(rows: usize, cols: usize, data: Vec<T>) -> Self
The matrix of rows × cols elements given row after row.
§Panics
Panics if data doesn’t hold rows · cols elements.
Sourcepub fn from_fn(
rows: usize,
cols: usize,
f: impl FnMut(usize, usize) -> T,
) -> Self
pub fn from_fn( rows: usize, cols: usize, f: impl FnMut(usize, usize) -> T, ) -> Self
The matrix whose element in row i and column j is f(i, j), computed column after
column.
§Panics
Panics if rows · cols overflows usize.
Sourcepub fn get(&self, i: usize, j: usize) -> Option<&T>
pub fn get(&self, i: usize, j: usize) -> Option<&T>
The element in row i and column j, or None outside the matrix.
Sourcepub fn get_mut(&mut self, i: usize, j: usize) -> Option<&mut T>
pub fn get_mut(&mut self, i: usize, j: usize) -> Option<&mut T>
The element in row i and column j, to change, or None outside the matrix.
Sourcepub fn column_mut(&mut self, j: usize) -> &mut [T]
pub fn column_mut(&mut self, j: usize) -> &mut [T]
Sourcepub fn row(&self, i: usize) -> impl ExactSizeIterator<Item = &T> + '_
pub fn row(&self, i: usize) -> impl ExactSizeIterator<Item = &T> + '_
Sourcepub fn as_mut_slice(&mut self) -> &mut [T]
pub fn as_mut_slice(&mut self) -> &mut [T]
The elements, column after column, to change.
Source§impl<T: Clone> Matrix<T>
impl<T: Clone> Matrix<T>
Trait Implementations§
Source§impl<T: Element> AddAssign for Matrix<T>
impl<T: Element> AddAssign for Matrix<T>
Source§fn add_assign(&mut self, rhs: Matrix<T>)
fn add_assign(&mut self, rhs: Matrix<T>)
+= operation. Read moreSource§impl<T: Element> AddAssign<&Matrix<T>> for Matrix<T>
impl<T: Element> AddAssign<&Matrix<T>> for Matrix<T>
Source§fn add_assign(&mut self, rhs: &Matrix<T>)
fn add_assign(&mut self, rhs: &Matrix<T>)
+= operation. Read moreSource§impl<T: Scalar> DivAssign<T> for Matrix<T>
impl<T: Scalar> DivAssign<T> for Matrix<T>
Source§fn div_assign(&mut self, rhs: T)
fn div_assign(&mut self, rhs: T)
/= operation. Read moreimpl<T: Eq> Eq for Matrix<T>
Source§impl<T: Element> MulAssign for Matrix<T>
impl<T: Element> MulAssign for Matrix<T>
Source§fn mul_assign(&mut self, rhs: Matrix<T>)
fn mul_assign(&mut self, rhs: Matrix<T>)
*= operation. Read moreSource§impl<T: Element> MulAssign<&Matrix<T>> for Matrix<T>
impl<T: Element> MulAssign<&Matrix<T>> for Matrix<T>
Source§fn mul_assign(&mut self, rhs: &Matrix<T>)
fn mul_assign(&mut self, rhs: &Matrix<T>)
*= operation. Read moreSource§impl<T: Element> MulAssign<&T> for Matrix<T>
impl<T: Element> MulAssign<&T> for Matrix<T>
Source§fn mul_assign(&mut self, rhs: &T)
fn mul_assign(&mut self, rhs: &T)
*= operation. Read moreSource§impl<T: Element> MulAssign<T> for Matrix<T>
impl<T: Element> MulAssign<T> for Matrix<T>
Source§fn mul_assign(&mut self, rhs: T)
fn mul_assign(&mut self, rhs: T)
*= operation. Read more