#### Overview of this book

Mastering SciPy
Credits
www.PacktPub.com
Preface
Free Chapter
Numerical Linear Algebra
Interpolation and Approximation
Differentiation and Integration
Nonlinear Equations and Optimization
Initial Value Problems for Ordinary Differential Equations
Computational Geometry
Descriptive Statistics
Inference and Data Analysis
Mathematical Imaging
Index

## 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

`norm(A,numpy.inf)`

Sum of absolute values of entries in each row. Pick the largest value.

`norm(A,-numpy.inf)`

Sum of absolute values of entries in each row. Pick the smallest value.

`norm(A,1)`

Sum of absolute values of entries in each column. Pick the largest value.

`norm(A,-1)`

Sum of absolute values of entries in each column. Pick the smallest value.

`norm(A,2)`

Largest eigenvalue of the matrix.

`norm(A,-2)`

Smallest eigenvalue of the matrix.

`norm(A,'fro')` or `norm(A,'f')`

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
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)