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
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 elementwise 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]])
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 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 nonsingular 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 nonsingular 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 MoorePenrose pseudoinverse 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 pseudoinverse. The function pinv2
computes the pseudoinverse 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 nonsingular matrices, the inverse and pseudoinverse 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
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 1norm 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 2norms, 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 singleprecision 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 nonsingular 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 wellconditioned). 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 illconditioned).
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 predefined 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