Basic matrix manipulation
The emphasis of the second part of this chapter is on mastering the following operations:
- Scalar multiplication, matrix addition, and matrix multiplication
- Traces and determinants
- Transposes and inverses
- Norms and condition numbers
Scalar multiplication, matrix addition, and matrix multiplication
Let us start with the matrices stored with the ndarray class. We accomplish scalar multiplication with the * operator, and the matrix addition with the + operator. But for matrix multiplication we will need the instance method dot() or the numpy.dot function, since the operator * is reserved for element-wise multiplication:
In [54]: 2*A Out[54]: array([[ 2, 4], [ 8, 32]]) In [55]: A + 2*A Out[55]: array([[ 3, 6], [12, 48]]) In [56]: A.dot(2*A) In [56]: np.dot(A, 2*A) Out[56]: Out[56]: array([[ 18, 68], array([[ 18, 68], [136, 528]]) [136, 528]]) In [57]: A.dot(B) ValueError: objects are not aligned In [58]: B.dot(A) In [58]: np.dot(B, A) Out[58]: Out[58]: array([[ -9, -34], array([[ -9, -34], [ 0, 0], [ 0, 0], [ 9, 34]]) [ 9, 34]])
The matrix class makes matrix multiplication more intuitive: the operator * can be used instead of the dot() method. Note also how matrix multiplication between different instance classes ndarray and a matrix is always casted to a matrix instance class:
In [59]: C * B ValueError: shapes (2,2) and (3,2) not aligned: 2 (dim 1) != 3 (dim 0) In [60]: B * C Out[60]: matrix([[ -9, -34], [ 0, 0], [ 9, 34]])
For sparse matrices, both scalar multiplication and addition work well with the obvious operators, even if the two sparse classes are not the same. Note the resulting class casting after each operation:
In [61]: S_10_coo = spsp.rand(5, 5, density=0.1, format='coo') In [62]: S_25_lil + S_10_coo Out[62]: <5x5 sparse matrix of type '<type 'numpy.float64'>' with 8 stored elements in Compressed Sparse Row format> In [63]: S_25_lil * S_10_coo Out[63]: <5x5 sparse matrix of type '<type 'numpy.float64'>' with 4 stored elements in Compressed Sparse Row format>
Tip
numpy.dot does not work well for matrix multiplication of a sparse matrix with a generic. We must use the operator * instead.
In [64]: S_100_coo = spsp.rand(2, 2, density=1, format='coo') In [65]: np.dot(A, S_100_coo) Out[66]: array([[ <2x2 sparse matrix of type '<type 'numpy.float64'>' with 4 stored elements in COOrdinate format>, <2x2 sparse matrix of type '<type 'numpy.float64'>' with 4 stored elements in COOrdinate format>], [ <2x2 sparse matrix of type '<type 'numpy.float64'>' with 4 stored elements in COOrdinate format>, <2x2 sparse matrix of type '<type 'numpy.float64'>' with 4 stored elements in COOrdinate format>]], dtype=object) In [67]: A * S_100_coo Out[68]: array([[ 1.81 , 1.555], [ 11.438, 11.105]])
Traces and determinants
The traces of a matrix are the sums of the elements on the diagonals (assuming always increasing indices in both dimensions). For generic matrices, we compute them with the instance method trace(), or with the function numpy.trace:
In [69]: A.trace() In [71]: C.trace() Out[69]: 17 Out[71]: matrix([[17]]) In [70]: B.trace() In [72]: np.trace(B, offset=-1) Out[70]: -1 Out[72]: 2
In order to compute the determinant of generic square matrices, we need the function det in the module scipy.linalg:
In [73]: spla.det(C) Out[73]: 8.0
Transposes and inverses
Transposes can be computed with any of the two instance methods transpose() or T, for any of the two classes of generic matrices:
In [74]: B.transpose() In [75]: C.T Out[74]: Out[75]: array([[-1, 0, 1], matrix([[ 1, 4], [-2, 0, 2]]) [ 2, 16]])
Hermitian transpose can be computed for the matrix class with the instance method H:
In [76]: D = C * np.diag((1j,4)); print D In [77]: print D.H [[ 0.+1.j 8.+0.j] [[ 0.-1.j 0.-4.j] [ 0.+4.j 64.+0.j]] [ 8.-0.j 64.-0.j]]
Inverses of non-singular square matrices are computed for the ndarray class with the function inv in the module scipy.linalg. For the matrix class, we may also use the instance method I. For non-singular square sparse matrices, we may use the function inv in the module scipy.sparse.linalg.
Tip
Inverses of sparse matrices are seldom sparse. For this reason, it is not recommended to perform this operation with the scipy.sparse.inv function. One possible way to go around this issue is to convert the matrix to generic with the todense() instance method, and use scipy.linear.inv instead.
But due to the difficulty of inverting large matrices, it is often beneficial to compute approximations to the inverse, instead. The function spilu in the module scipy.sparse.linalg provides us with a very fast algorithm to perform this computation for square sparse matrices in CSC format. This algorithm is based on LU decompositions, and coded internally as a wrapper of a function from the library SuperLU. Its use is rather complex, and we are going to postpone its study until we explore matrix factorizations.
In [78]: E = spsp.rand(512, 512, density=1).todense() In [79]: S_100_csc = spsp.rand(512, 512, density=1, format='csc') In [80]: %timeit E.I 10 loops, best of 3: 28.7 ms per loop In [81]: %timeit spspla.inv(S_100_csc) 1 loops, best of 3: 1.99 s per loop
Note
In the execution of sparse inverses, if the input matrix is not in the CSC or CSR format, we will get a warning:
/scipy/sparse/linalg/dsolve/linsolve.py:88: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format warn('spsolve requires A be CSC or CSR matrix format', SparseEfficiencyWarning) /scipy/sparse/linalg/dsolve/linsolve.py:103: SparseEfficiencyWarning: solve requires b be CSC or CSR matrix format
The Moore-Penrose pseudo-inverse can be computed for any kind of matrix (not necessarily square) with either routines the pinv or the pinv2 in the module scipy.linalg. The first method, pinv, resorts to solving a least squares problem to compute the pseudo-inverse. The function pinv2 computes the pseudo-inverse by a method based on singular value decompositions. For Hermitian matrices, or matrices that are symmetric with no complex coefficients, we also have a third function called pinvh, which is based on eigenvalue decompositions.
It is known that in the case of square non-singular matrices, the inverse and pseudo-inverse are the same. This simple example shows the times of computation of the inverses of a large generic symmetric matrix with the five methods described:
In [82]: F = E + E.T # F is symmetric In [83]: %timeit F.I 1 loops, best of 3: 24 ms per loop In [84]: %timeit spla.inv(F) 10 loops, best of 3: 28 ms per loop In [85]: %timeit spla.pinvh(E) 1 loops, best of 3: 120 ms per loop In [86]: %timeit spla.pinv2(E) 1 loops, best of 3: 252 ms per loop In [87]: %timeit spla.pinv(F) 1 loops, best of 3: 2.21 s per loop
Norms and condition numbers
For generic matrices, we have seven different standard norms in scipy.linalg. We can summarize them in the following table:
|
Constructor |
Description |
|---|---|
|
|
Sum of absolute values of entries in each row. Pick the largest value. |
|
|
Sum of absolute values of entries in each row. Pick the smallest value. |
|
|
Sum of absolute values of entries in each column. Pick the largest value. |
|
|
Sum of absolute values of entries in each column. Pick the smallest value. |
|
|
Largest eigenvalue of the matrix. |
|
|
Smallest eigenvalue of the matrix. |
|
|
Frobenius norm: the square root of the trace of the product A.H * A. |
In [88]: [spla.norm(A,s) for s in (np.inf,-np.inf,-1,1,-2,2,'fro')] Out[88]: [20, 3, 5, 18, 0.48087417361008861, 16.636368595013604, 16.643316977093239]
Tip
For sparse matrices, we can always compute norms by applying the todense() instance method prior to computation. But when the sizes of the matrices are too large, this is very impractical. In those cases, the best we can get for the 1-norm is a lower bound, thanks to the function onenormest in the module scipy.sparse.linalg:
In [89]: spla.norm(S_100_csc.todense(), 1) - \ ....: spspla.onenormest(S_100_csc) Out[89]: 0.0
As for the 2-norms, we may find the values of the smallest and the largest eigenvalue, but only for square matrices. We have two algorithms in the module scipy.sparse.linalg that perform this task: eigs (for generic square matrices) and eigsh for real symmetric matrices. We will explore them in detail when we discuss matrix decompositions and factorizations in the next section.
Note the subtle difference between the norm computations from SciPy and NumPy. For example, in the case of the Frobenius norm, scipy.linalg.norm is based directly on the BLAS function called NRM2, while numpy.linalg.norm is equivalent to a purely straightforward computation of the form sqrt(add.reduce((x.conj() * x).real)). The advantage of the code based on BLAS, besides being much faster, is clear when some of the data is too large or too small in single-precision arithmetic. This is shown in the following example:
In [89]: a = np.float64([1e20]); \ ....: b = np.float32([1e20]) In [90]: [np.linalg.norm(a), spla.norm(a)] Out[90]: [1e+20, 1e+20] In [91]: np.linalg.norm(b) [...]/numpy/linalg/linalg.py:2056: RuntimeWarning: overflow encountered in multiply return sqrt(add.reduce((x.conj() * x).real, axis=None)) Out[91]: inf In [92]: spla.norm(b) Out[92]: 1.0000000200408773e+20
This brings us inevitably to a discussion about the computation of the condition number of a non-singular square matrix A. This value measures how much the output of the solution to the linear equation A * x = b will change when we make small changes to the input argument b. If this value is close to one, we can rest assured that the solution is going to change very little (we say then that the system is well-conditioned). If the condition number is large, we know that there might be issues with the computed solutions of the system (and we say then that it is ill-conditioned).
The computation of this condition number is performed by multiplying the norm of A with the norm of its inverse. Note that there are different condition numbers, depending on the norm that we choose for the computation. These values can also be computed for each of the pre-defined norms with the function numpy.linalg.cond, although we need to be aware of its obvious limitations.
In [93]: np.linalg.cond(C, -np.inf) Out[93]: 1.875
