In [2]:

```
from numpy import identity, savetxt
from scipy.io import savemat
n = 1000
denseI = identity(n, dtype='float64')
savetxt('tmp/denseI.txt', denseI)
dense_dict = {'denseI': denseI}
savemat('tmp/denseI.mat', dense_dict)
```

In [3]:

```
!ls -l tmp/denseI*
```

In [4]:

```
from scipy.sparse import identity
sparseI = identity(n, dtype='float64', format='csr')
sparse_dict = {'sparseI' : sparseI}
savemat('tmp/sparseI.mat', sparse_dict)
```

In [5]:

```
!ls -l tmp/sparseI*
```

In [6]:

```
from numpy import array
from scipy.sparse import coo_matrix
rows = array([0, 0, 1, 2, 2, 2])
cols = array([0, 2, 2, 0, 1, 2])
data = array([1, 2, 3, 4, 5, 6])
C = coo_matrix((data, (rows, cols)), shape=(3,3))
C.todense() # converting to NumPy matrix to print contents easily
```

Out[6]:

In [7]:

```
C
```

Out[7]:

In [8]:

```
C.nnz
```

Out[8]:

In [9]:

```
C.T.todense()
```

Out[9]:

In [10]:

```
from scipy.sparse import csr_matrix
indptr = array([0,2,3,6]) # row 0 indices start at index 0 of indices,
# row 1 indices start at index 2 of indices
# row 3 indices start at index 3 of indices
# Last item is length of data
indices = array([0,2,2,0,1,2]) # column numbers
data = array([1,2,3,4,5,6])
R = csr_matrix((data, indices, indptr), shape=(3,3))
R.todense()
```

Out[10]:

In [11]:

```
R.tocoo()
```

Out[11]:

In [12]:

```
R
```

Out[12]:

In [13]:

```
C.tocsr()
```

Out[13]:

In [14]:

```
C
```

Out[14]:

In [15]:

```
from scipy.sparse import lil_matrix
L = lil_matrix((3,3), dtype='float64')
```

In [16]:

```
L
```

Out[16]:

In [17]:

```
L.todense()
```

Out[17]:

In [18]:

```
L
```

Out[18]:

In [19]:

```
L[0,0] = 2; L[1,2] = 3
```

In [20]:

```
L.todense()
```

Out[20]:

In [21]:

```
L[1,1:3].todense()
```

Out[21]:

In [22]:

```
L[1,1:3]
```

Out[22]:

In [23]:

```
C[1,1:3]
```

In [24]:

```
R[1,1:3]
```

Out[24]:

In [25]:

```
C[0,0] = 2; C[1,2] = 3
```

In [26]:

```
R[0,0] = 111; R[1,2] = 999
R.todense()
```

Out[26]:

- Write the CSR representation for the matrix: $$ A= \left[ \begin{array}{ccc} 1 & 0 & 3 \\ 0 & 5 & 0 \\ 0 & 0 & 2 \end{array} \right] $$ and then use that to create a CSR matrix. Then create the same matrix by first creating the appropriate COO matrix.
- Construct a dense 5 x 5 matrix of floats with 2 along the main diagonal, and -1 along the diagonal above and below the main diagonal. We will refer to this type of matrix as a discrete 1D Laplacian (we'll see the reason for the name in the differential equation unit). You can use
`numpy.eye`

or`scipy.linalg.toeplitz`

. If you use the`eye`

then you can build identity matrices using`eye(n)`

and zero matrices with an off-diagonal of ones using`eye(n, k=...|-1|0|1|...)`

.

`scipy.linalg`

without worrying about the choices available amongst the solving techniques. We will cover that below. First suppose you have a small dense `numpy`

array and just want to solve a linear system once.

In [2]:

```
from scipy.linalg import solve
import numpy as np
A = np.array(range(1,9) + [0], dtype=float).reshape((3,3))
print A
b = np.array([1., 0, 2])
x = solve(A, b)
print
print x
```

*condition number* of a matrices that we will look at later. But for now lets just compute $A\; x$ using the solution $x$ obtained above.

In [3]:

```
np.dot(A, x)
```

Out[3]:

`scipy.linalg`

module has a function `lu`

which takes an array and returns its LU decomposition. In addition, the module provides a function `solve_triangular(T, b, lower=True|False)`

which efficiently solves triangular matrices. We can combine these to solve the linear system we solved abovebut this time using the procedure outlined above.

In [27]:

```
import numpy as np
```

In [5]:

```
from scipy.linalg import lu, solve_triangular
A = np.array(range(1,9) + [0], dtype=float).reshape((3,3))
b = np.array([1., 0, 2])
P, L, U = lu(A)
y = solve_triangular(L, np.dot(P.T, b), lower=True)
x = solve_triangular(U, y)
print x
```

`scipy.linalg`

for matrices with other commonly occuring structure. Many matrices are banded, in that there is a band of nonzero entries only in some banded neigborhood of the diagonal. Correspondingly there is `solve_banded`

function for this. For Hermitian (or symmetric in the real case) matrices there are special functions. You should explore the `scipy.linalg`

reference pages.

Lets take another matrix and examine the factors $P$, $L$, and $U$.

In [28]:

```
A = np.array([[2,1,1,0], [4,3,3,1], [8,7,9,5], [6,7,9,8]], dtype=float)
print A
```

In [29]:

```
from scipy.linalg import lu
P, L, U = lu(A)
print L
```

In [30]:

```
print U
```

In [31]:

```
print P
```

$P$ is a permutation of the identity matrix of order 4 and $P^T\, P = P\, P^T = I$.

In [44]:

```
np.dot(P.T, P)
```

Out[44]:

In [32]:

```
print np.dot(L, U)
```

In [33]:

```
np.dot(P, np.dot(L, U))
```

Out[33]:

In [34]:

```
print A
```

*pivoting* which involves the use of the largest entry in each column and swapping rows.

*Cholesky factorization* of $A$:
$$
A = L \, L^T\, .
$$
Since pivoting is not needed, there is no permutation matrix produced in the case of Cholesky. Some people use Cholesky to discover if a symmetric matrix is positive definite.

In [35]:

```
from scipy.linalg import cholesky
A = np.array([[2.,-1,0], [-1,2,-1], [0,-1,2]])
print A
L = cholesky(A, lower=True)
np.dot(L, L.T)
```

Out[35]:

- Create three random 5 x 5 matrices $A, B$ and $C$ and pick a random vector $b$ with 5 entries. Compute $x = B^{-1}\, (2A + I)\,(C^{-1} + A)\; b$ without computing any matrix inverses.
- Use LU factorization solve a 100 x 100 random linear system using the procedure outlined above. To check the quality of the solution, first start with a known random vector
`xtrue`

and compute the RHS vector`b`

as`dot(A, xtrue)`

. This way we will have the true solution available to verify the answer. Now pretend you don't know the solution and solve $A\; x = b$ for unknown $x$ using the $b$ you generated above. Compare the resulting solution`x`

with the true solution`xtrue`

by computing the norm of the difference. This can be done using`numpy.linalg.norm`

.

*Jacobi iteration* and if $M$ is the lower triangular part of $A$ the iteration is called *Gauss-Seidel iteration*.

*conjugate gradient* method that should be the first thing to try for symmetric positive definite matrices and the *GMRES* method for general matrices. The details of these methods are beyond the scope of this but you should definitely learn to use these as available in `scipy.sparse.linalg`

.

- Implement Jacobi iteration for the 5 x 5 discrete 1D Laplacian matrix $A$ and solve the linear system with RHS being a vector of all ones. Run the interation for 50 steps and print out the norm of the residual as you go. Or you can implement reasonable stopping criteria such as stopping when progress stalls or some maximum number of iterations has been exceeded.
- Use the conjugate gradient solver
`cg`

from`scipy.sparse.linalg`

to solve the system $A\, x = b$ where $b$ is a vector of ones and $A$ is a sparse discrete 1D Laplacian matrix from the previous problem. Print out the residual at each step using the optional argument`callback=`

. This argument takes a function`f(xk)`

which is called during each step and provides the user with the current guess for the solution`xk`

.

The *condition number* of a square matrix $A$ is $\operatorname{cond}(A):= || A ||\; || A^{-1} ||$.

If instead of the RHS $b$ one perturbs the matrix $A$ and $E$ is this perturbation then

$$ \frac{||\Delta x||}{||\hat{x}||} \le \operatorname{cond}(A)\frac{||E||}{||A||} $$*residual*. The quantity $||r|| / (||A||\, ||\hat{x}||)$ is called the *relative residual*. The error and residual are related by
$$
\frac{||\Delta x||}{||\hat{x}||} \le \operatorname{cond}(A)\frac{||r||}{||A||\;||\hat{x}||}
$$

The *Hilbert matrix* is parametrized by its size and gets increasingly ill-conditioned with size.

In [7]:

```
from scipy.linalg import hilbert
print hilbert(4)
```

In [8]:

```
from numpy.linalg import cond
for n in [2, 5, 8, 11, 14]:
print cond(hilbert(n))
```

In [9]:

```
from numpy.random import rand, seed
seed(6172014)
for n in [2, 5, 8, 11, 14]:
print cond(rand(n, n))
```

The impact of condition number on the error is visible very dramatically for Hilbert matrices.

In [10]:

```
from numpy.linalg import solve, norm
import numpy as np
for n in [2, 5, 8, 11, 14]:
A = hilbert(n)
b = np.dot(A, np.ones(n))
x = solve(A, b)
abserror = norm(np.ones(n) - x)
relerror = abserror/norm(np.ones(n))
print '%2d\t%2.5e\t%2.5e\t%2.5e' % (n, abserror, relerror, norm(b-np.dot(A,x)))
```

For the random matrices, due to the smaller condition numbers the errors don't grow so badly.

In [11]:

```
seed(6172014)
for n in [2, 5, 8, 11, 14]:
A = rand(n, n)
b = np.dot(A, np.ones(n))
x = solve(A, b)
abserror = norm(np.ones(n) - x)
relerror = abserror/norm(np.ones(n))
print '%2d\t%2.5e\t%2.5e\t%2.5e' % (n, abserror, relerror, norm(b-np.dot(A,x)))
```

- A square matrix is called
*orthogonal*if it has orthogonal columns. How would you express that as a matrix equation? What is the 2-norm of an orthogonal matrix? What is its condition number?

To be continued ...