Linear Algebra

`MyMatrix`

A matrix utility class that provides a set of methods for matrix operations.

Attributes:

Name	Type	Description
`X (np.ndarray)`	`A matrix.`
`row (int)`	`Number of rows.`
`col (int)`	`Number of columns.`
`det (np.ndarray)`	`The determinant of X, if it is square.`

`init(X, axis=None)`

Initialize the MyMatrix instance.

Parameters:

Name	Type	Description	Default
`X`	`Union[ndarray, List[List[Any]], Tuple[Tuple[Any]]]`	A matrix or a sequence that represents the matrix.	required
`axis`	`Optional[int]`	The axis along which the sequence is stacked into a matrix. Required if X is a sequence.	`None`

Returns:

Type	Description
`None`

`repr()`

Delegate to numpy's representation.

`approx_rank_k(k)`

Computes the rank-k approximation of the matrix X. The matrix approximation can be constructed as the matrix product u @ np.diag(s) @ vh where (u, s, vh) are the tuple returned by calling this method.

Parameters:

Name	Type	Description	Default
`k`	`int`	The rank of the approximation.	required

Returns:

Type	Description
`Tuple[ndarray, ndarray, ndarray]`	The singular values, left, right singular vectors needed to construct the rank-k approximation of the matrix X.

Examples:

>>> import linalg as la
>>> m = la.MyMatrix([[1, 0, 12], [7, 12, 17], [9, 77, 27], [8, 7, 16]], axis=0)
>>> m.approx_rank_k(1)
(array([[-0.05593827, -0.49093673],
        [-0.2166436 , -0.55228019],
        [-0.96140112,  0.25677546],
        [-0.16013852, -0.62292382]]),
 array([85.19359971, 21.49336595]),
 array([[-0.13505899, -0.91261056, -0.38587696],
        [-0.32704554,  0.40867874, -0.85206978]]))

`char_poly()`

Computes characteristic polynomial det(lambda*I - M) where I is the identity matrix.

Returns:

Type	Description
`sympy expression`	A class instance inheriting from the Base class for algebraic expressions `Expr`.

Raises:

Type	Description
`LinAlgError`	The matrix X must be square.

Examples:

>>> import linalg as la
>>> m = la.MyMatrix([[1, 2], [3, 4]], axis=0)
>>> m.char_poly()
lambda**2 - 5*lambda - 2

`cob_mat(original_basis, new_basis, inv=False)` `staticmethod`

Find the change of basis matrix S from original_basis to new_basis. The change of basis matrix from basis A to basis B is defined to be: S_{A -> B} = [ [a_1]_B, ..., [a_n]_B ].

Parameters:

Name	Type	Description	Default
`original_basis`	`List[List]`	A set of `n` basis vectors for :math:`R^{n}` or `m` basis vectors for a subspace of :math:`R^{n}`.	required
`new_basis`	`List[List]`	A set of `n` basis vectors for :math:`R^{n}` or `m` basis vectors for a subspace of :math:`R^{n}`.	required
`inv`	`bool`	Whether to return the inverse of the change of basis matrix, by default False.	`False`

Returns:

Type	Description
`ndarray`	The change of basis matrix or the inverse of it, which translates from `new_basis` back to `original_basis`.

Examples:

>>> import linalg as la
>>> # Change of basis from [[1, 2, -3], [4, -1, -3]] to [[0, 1, -1], [1, 0, -1]]
>>> la.MyMatrix.cob_mat([[1, 2, -3], [4, -1, -3]], [[0, 1, -1], [1, 0, -1]], False)
array([[ 2., -1.],
       [ 1.,  4.]])

`col_space()`

Returns a list of vectors (np.ndarray objects) that span the column-space or image of X.

Returns:

Type	Description
`List[ndarray]`	A list of columns vectors.

Examples:

>>> import linalg as la
>>> m = la.MyMatrix([[1, 2], [3, 4], [5, 6]], axis=0)
>>> [array([[1.],
            [3.],
            [5.]]),
    array([[2.],
            [4.],
            [6.]])]

`diagonalize(reals_only=False, sort=False, normalize=False)`

Return (P, D), where D is a diagonal D = P^-1 * M * P where M is current matrix.

Parameters:

Name	Type	Description	Default
`reals_only`	`bool`	Whether to throw an error if complex numbers are need to diagonalize, by default False.	`False`
`sort`	`bool`	Whether to sort the eigenvalues along the diagonal, by default False.	`False`
`normalize`	`bool`	Whether to normalize the eigenvector columns of P, by default False.	`False`

Returns:

Type	Description
`Tuple[ndarray, ndarray]`	A tuple (P, D).

Raises:

Type	Description
`LinAlgError`	The matrix X must be square.

Examples:

>>> import linalg as la
>>> m = la.MyMatrix([[1, 2], [2, 1]], axis=1)
>>> m.diagonalize()
(array([[-1.,  1.],
        [ 1.,  1.]]),
 array([[-1.,  0.],
        [ 0.,  3.]]))

`eigen(real=True)`

Find the eigenvalues and eigenvectors of the matrix X.

Parameters:

Name	Type	Description	Default
`real`	`bool`	Whether to convert the array of eigenvalues to real numbers, by default True.	`True`

Returns:

Type	Description
`w (…, M) array`	The eigenvalues, each repeated according to its multiplicity. The eigenvalues are not necessarily ordered. The resulting array will be of complex type, unless the imaginary part is zero in which case it will be cast to a real type. When `X` is real the resulting eigenvalues will be real (0 imaginary part) or occur in conjugate pairs.
`vl(…, M, M) array`	The normalized (unit “length”) left eigenvectors, such that the column v[:,i] is the eigenvector corresponding to the eigenvalue w[i].
`vr(…, M, M) array`	The normalized (unit “length”) right eigenvectors, such that the column v[:,i] is the eigenvector corresponding to the eigenvalue w[i].

Raises:

Type	Description
`LinAlgError`	The matrix X must be square.

Examples:

>>> import linalg as la
>>> m = la.MyMatrix([[1, 2], [2, 1]], axis=0)
>>> m.eigen()
(array([ 3., -1.]),
 array([[ 0.70710678, -0.70710678],
        [ 0.70710678,  0.70710678]]),
 array([[ 0.70710678, -0.70710678],
        [ 0.70710678,  0.70710678]]))

`eye(N, M=None, k=0, dtype=np.dtype(np.int8))` `classmethod`

Instantiate a class instance with an identity matrix.

Parameters:

Name	Type	Description	Default
`N`	`int`	Number of rows.	required
`M`	`Optional[int]`	Number of columns, by default None, which takes the value of N.	`None`
`k`	`int`	Index of the diagonal where 0 (the default) refers to the main diagonal, a positive value refers to an upper diagonal, and a negative value to a lower diagonal, by default 0.	`0`
`dtype`	`Optional[dtype]`	Data-type of identity matrix, by default np.int8.	`dtype(int8)`

Returns:

Type	Description
`MyMatrix`	An instance of of class MyMatrix.

Examples:

>>> import linalg as la
>>> la.MyMatrix.eye(3)
array([[1, 0, 0],
       [0, 1, 0],
       [0, 0, 1]], dtype=int8)
>>> la.MyMatrix.eye(3, 4, 1)
array([[0, 1, 0, 0],
       [0, 0, 1, 0],
       [0, 0, 0, 1]], dtype=int8)

`gs(ret_type='matrix')`

Create an orthogonal matrix (a set of orthonormal basis vectors) for matrix X.

Parameters:

Name	Type	Description	Default
`ret_type`	`str`	Whether to return a set of orthonormal basis vectors or an orthogonal matrix, by default 'matrix'.	`'matrix'`

Returns:

Type	Description
`ndarray or MyMatrix`	A set of orthonormal basis vectors or an orthogonal matrix, depending on the `ret_type` argument.

Raises:

Type	Description
`LinAlgError`	The matrix X must be square.

Examples:

>>> import linalg as la
>>> m = la.MyMatrix([[1, 0], [0, 1], [1, 1]], axis=0)
>>> m.gs(ret_type="vector")
[array([-0.70710678, -0.        , -0.70710678]), array([ 0.40824829, -0.81649658, -0.40824829])]
>>> m.gs(ret_type="matrix")
array([[-0.70710678,  0.40824829],
       [-0.        , -0.81649658],
       [-0.70710678, -0.40824829]])

`inv()`

Given a square matrix X, return the matrix Xinv satisfying dot(X, Xinv) = dot(Xinv, X) = eye(X.shape[0]).

Returns:

Type	Description
`MyMatrix`	The inverse of X.

Raises:

Type	Description
`LinAlgError`	If `X` is not square or inversion fails.

Examples:

>>> import linalg as la
>>> m = la.MyMatrix([[1, 2], [3, 4]], axis=0)
>>> m.inv()
>>> [array([-2.,  1.])]
           [ 1.5, -0.5]])

`is_coord_vec(original_vec, coord_vec, *args)` `staticmethod`

Test if coord_vec is the coordinate vector with respect to a set of basis vectors basis_vecs.

Parameters:

Name	Type	Description	Default
`original_vec`	`Union[ndarray, List, Tuple]`	The original vector in its original coordinate.	required
`coord_vec`	`Union[ndarray, List, Tuple]`	A coordinate vector.	required
`args`	`Union[ndarray, List, Tuple]`	An arbitrary number of basis vectors.	`()`

Returns:

Type	Description
`bool`	Whether `coord_vec` is the coordinate vector with respect to the basis vectors.

Examples:

>>> import linalg as la
>>> original_vec = np.array([5, 11])
>>> coord_vec = [1, 2]
>>> basis_vecs = [np.array([1, 3]), np.array([2, 4])]
>>> la.MyMatrix.is_coord_vec(original_vec, coord_vec, *basis_vecs)
True

`is_pos_def()`

Check if matrix X is positive definite. A symmetric matrix A is positive definite if (and only if) all of its eigenvalues are positive. The matrix A is positive sem-definite if (and only if) all of its eigenvalues are non-negative (positive or zero).

Returns:

Type	Description
`bool`	True if positive definite or False if non-positive definite.

Examples:

>>> import linalg as la
>>> m = la.MyMatrix.eye(3)
>>> m.is_pos_def()
True

`is_sq()`

Check if matrix X is square.

Returns:

Type	Description
`bool`	True if square or False if non-square.

Examples:

>>> import linalg as la
>>> la.MyMatrix(np.array([[1, 2], [3, 4]]), axis=0).is_sq()
True
>>> la.MyMatrix(np.array([[1, 2], [3, 4], [5, 6]]), axis=0).is_sq()
False

`null_space()`

Returns a list of vectors (np.ndarray objects) that span the null-space or kernel of X.

Returns:

Type	Description
`List[ndarray]`	A list of columns vectors.

Examples:

>>> import linalg as la
>>> m = la.MyMatrix([[1, 2], [3, 6]], axis=0)
>>> m.null_space()
>>> [array([[-2.],
            [ 1.]])]

`pinv()`

Compute the (Moore-Penrose) pseudo-inverse of a matrix. Calculate the generalized inverse of a matrix using its singular-value decomposition (SVD) and including all large singular values.

Returns:

Name	Type	Description
`B`	`(..., N, M) MyMatrix`	The pseudo-inverse of `X`.

Raises:

Type	Description
`LinAlgError`	If the SVD computation does not converge.

Examples:

>>> import linalg as la
>>> m = la.MyMatrix([[1, 2], [3, 4], [5, 6]], axis=1)
>>> m.pinv()
array([[-1.33333333,  1.08333333],
    [-0.33333333,  0.33333333],
    [ 0.66666667, -0.41666667]])

`power(n)`

Raise a square matrix to the (integer) power n.

Parameters:

Name	Type	Description	Default
`n`	`int`	The exponent can be any integer or long integer, positive, negative, or zero.	required

Returns:

Type	Description
`MyMatrix`	An instance of class MyMatrix.

Raises:

Type	Description
`LinAlgError`	For matrices that are not square or that (for negative powers) cannot be inverted numerically.

Examples:

>>> import linalg as la
>>> la.MyMatrix(np.array([[1, 2], [12, 17]]), axis=0).power(2)
array([[ 25,  36],
       [216, 313]])

`rapprox_rank_k(k, n_oversamples=None, n_iters=None, return_onb=True)`

Computes the rank-k approximation of the matrix X using randomized SVD. The matrix approximation can be constructed as the matrix product u @ np.diag(s) @ vh where (u, s, vh) are the tuple returned by calling this method.

Parameters:

Name	Type	Description	Default
`k`	`int`	The rank of the approximation.	required
`n_oversamples`	`Optional[int]`	Additional number of random vectors to sample the column space of X so as to ensure proper conditioning.	`None`
`n_iters`	`Optional[int]`	Number of power iterations.	`None`
`return_onb`	`Optional[bool]`	Whether to return the orthonormal basis Q for the approximated column space of X.	`True`

Returns:

Type	Description
`Union[Tuple[ndarray, ndarray, ndarray], Tuple[ndarray, ndarray, ndarray, ndarray]]`	The singular values, left, right singular vectors needed to construct the rank-k approximation of the matrix X. Optionally, the orthonormal basis Q for the approximated column space of X can be returned.

Examples:

>>> import linalg as la
>>> m = la.MyMatrix([[1, 0, 12], [7, 12, 17], [9, 77, 27], [8, 7, 16]], axis=0)
>>> m.rapprox_rank_k(2)
(array([[-0.05593827,  0.49093673],
        [-0.2166436 ,  0.55228019],
        [-0.96140112, -0.25677546],
        [-0.16013852,  0.62292382]]),
 array([85.19359971, 21.49336595]),
 array([[-0.13505899, -0.91261056, -0.38587696],
        [ 0.32704554, -0.40867874,  0.85206978]]),
 array([[-0.58506352, -0.35851766, -0.70637316,  0.17378931],
        [-0.31208547, -0.44629115,  0.29152198, -0.78641071],
        [ 0.67509565, -0.71187709, -0.18086095,  0.06903766],
        [-0.3233407 , -0.40684188,  0.61914555,  0.58871833]]))

`rref(pivots=False)`

Return the reduced row-echelon form of the matrix. Use pivots to return the indices of pivot columns.

Parameters:

Name	Type	Description	Default
`pivots`	`bool`	Whether to return a tuple of pivot indices, by default False.	`False`

Returns:

Type	Description
`MyMatrix or tuple`	The reduced row-echelon from or a tuple of pivot column indices.

Examples:

>>> import linalg as la
>>> m = la.MyMatrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]], axis=1)
>>> m.rref()
array([[1., 0., -1.],
       [0., 1.,  2.],
       [0., 0.,  0.]])
>>> m.rref(pivots=True)
(array([[1., 0., -1.],
        [0., 1.,  2.],
        [0., 0.,  0.]]), (0, 1))

`sv_plot()`

Plot the log and cumulative sum of singular values of the matrix X. This can be used to visually assess how much information is captured by the first k-rank of the matrix X, so that a sensible number can be selected for k for rank-k matrix approximation.

Returns:

Type	Description
`None`

Examples:

>>> import linalg as la
>>> m = la.MyMatrix([[1, 0, 12], [7, 12, 17], [9, 77, 27], [8, 7, 16]], axis=0)
>>> m.sv_plot()

`svd(full_matrices=True, sigma_only=False)`

Factorizes the matrix X into two unitary matrices U and Vh, and a 1-D array s of singular values (real, non-negative) such that a == U @ S @ Vh, where S is a suitably shaped matrix of zeros with main diagonal s.

Parameters:

Name	Type	Description	Default
`full_matrices`	`Optional[bool]`	If True, u and vh have the shapes (..., M, M) and (..., N, N), respectively. Otherwise, the shapes are (..., M, K) and (..., K, N), respectively, where K = min(M, N)., by default True.	`True`
`sigma_only`	`Optional[bool]`	Whether to return the singular values only, by default False, which constructs the sigma matrix in SVD from the singular values.	`False`

Returns:

Type	Description
`Tuple[ndarray, ndarray, ndarray]`	A tuple of three matrices: U (m x m), S (m x n), V^T (n x n).

Examples:

>>> import linalg as la
>>> m = la.MyMatrix([[1, 0], [71, 1], [12, 1]], axis=1)
>>> m.svd()
(array([[-0.99987192, -0.01600465],
        [-0.01600465,  0.99987192]]),
 array([[72.02311126+0.j,  0.        +0.j,  0.        +0.j],
        [ 0.        +0.j,  0.81941679+0.j,  0.        +0.j]]),
 array([[-0.01388265, -0.98589063, -0.16681406],
        [-0.01953176, -0.16653093,  0.98584277],
        [-0.99971285,  0.01694429, -0.01694429]]))

`zeros(shape, dtype=np.dtype(np.float64))` `classmethod`

Instantiate a class instance with a zero matrix.

Parameters:

Name	Type	Description	Default
`shape`	`Tuple[int]`	Shape of the zero matrix, e.g. (2, 4).	required
`dtype`	`Optional[dtype]`	The desired data-type for the matrix, by default np.float64.	`dtype(float64)`

Returns:

Type	Description
`MyMatrix`	An instance of of class MyMatrix.

Examples:

>>> import linalg as la
>>> la.MyMatrix.zeros((2, 3))
array([[0., 0., 0.],
       [0., 0., 0.]])

`cob(vector, *args)`

Express a vector in a new basis spanned by a pairwise orthogonal (not orthonomal) set basis vectors.

Parameters:

Name	Type	Description	Default
`vector`	`Union[ndarray, List, Tuple]`	An n-dimensional vector.	required
`args`	`Union[ndarray, List, Tuple]`	An arbitrary number of pairwise orthogonal n-dimensional basis vectors.	`()`

Returns:

Type	Description
`ndarray`	The original vector expressed in a new basis.

Examples:

>>> import linalg as la
>>> # Express [10, -5] in the basis spanned by [3, 4] and [4, -3]
>>> la.cob(vector=np.array([10, -5]), np.array([3, 4]), np.array([4, -3]))

`lin_ind(*args)`

Check if the input vectors are linearly independent. If not, return the redundant column vectors as an np.ndarray.

Parameters:

Name	Type	Description	Default
`args`	`Union[ndarray, List, Tuple]`	An arbitrary number of vectors.	`()`

Returns:

Type	Description
`None or ndarray`	An `np.ndarray` of redundant columns or None if the input vectors are linearly independent.

Examples:

>>> import linalg as la
>>> # Linear independent
>>> la.lin_ind([1, 0, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1], [0, 1, 0, 0])
The matrix with input vectors as columns has full column rank, and so the input vectors are linearly independent
>>> # Linear dependent
>>> lin_ind([1, 0, 0, 0], [1, 9, 3, 0], [0, 0, 4, 1], [0, 1, 0, 0], [2, 3, 4, 9])
array([[2, 3, 4, 9]])

`proj(v1, v2, proj_type='vector')`

Obtain the scalar or vector projection of vector v1 onto vector v2.

Parameters:

Name	Type	Description	Default
`v1`	`Union[ndarray, List, Tuple]`	An n-dimensional vector.	required
`v2`	`Union[ndarray, List, Tuple]`	An n-dimensional vector.	required
`proj_type`	`str`	Type of projection, by default 'vector'.	`'vector'`

Returns:

Type	Description
`ndarray`	The scalar or vector projection of vector `v1` onto vector `v2`.

Examples:

>>> import linalg as la
>>> la.proj(np.array([2, 4, 0]), np.array([4, 2, 4]), 'scalar')
np.float64(2.6666666666666665)
>>> la.proj(np.array([2, 1]), np.array([3, -4]), 'vector')
array([ 0.24, -0.32])

Linear Algebra

MyMatrix

__init__(X, axis=None)

__repr__()

approx_rank_k(k)

char_poly()

cob_mat(original_basis, new_basis, inv=False) staticmethod

col_space()

diagonalize(reals_only=False, sort=False, normalize=False)

eigen(real=True)

eye(N, M=None, k=0, dtype=np.dtype(np.int8)) classmethod

gs(ret_type='matrix')

inv()

is_coord_vec(original_vec, coord_vec, *args) staticmethod

is_pos_def()

is_sq()

null_space()

pinv()

power(n)

rapprox_rank_k(k, n_oversamples=None, n_iters=None, return_onb=True)

rref(pivots=False)

sv_plot()

svd(full_matrices=True, sigma_only=False)

zeros(shape, dtype=np.dtype(np.float64)) classmethod

cob(vector, *args)

lin_ind(*args)

proj(v1, v2, proj_type='vector')

`MyMatrix`

`init(X, axis=None)`

`repr()`

`approx_rank_k(k)`

`char_poly()`

`cob_mat(original_basis, new_basis, inv=False)` `staticmethod`

`col_space()`

`diagonalize(reals_only=False, sort=False, normalize=False)`

`eigen(real=True)`

`eye(N, M=None, k=0, dtype=np.dtype(np.int8))` `classmethod`

`gs(ret_type='matrix')`

`inv()`

`is_coord_vec(original_vec, coord_vec, *args)` `staticmethod`

`is_pos_def()`

`is_sq()`

`null_space()`

`pinv()`

`power(n)`

`rapprox_rank_k(k, n_oversamples=None, n_iters=None, return_onb=True)`

`rref(pivots=False)`

`sv_plot()`

`svd(full_matrices=True, sigma_only=False)`

`zeros(shape, dtype=np.dtype(np.float64))` `classmethod`

`cob(vector, *args)`

`lin_ind(*args)`

`proj(v1, v2, proj_type='vector')`