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

## Creation of matrices and linear operators

In the first part of this chapter, we are going to focus on the effective creation of matrices. We start by recalling some different ways to construct a basic matrix as an `ndarray` instance class, including an enumeration of all the special matrices already included in NumPy and SciPy. We proceed to examine the possibilities of constructing complex matrices from basic ones. We review the same concepts within the `matrix` instance class. Next, we explore in detail the different ways to input sparse matrices. We finish the section with the construction of linear operators.

### Note

We assume familiarity with `ndarray` creation in NumPy, as well as data types (dtype), indexing, routines for the combination of two or more arrays, array manipulation, or extracting information from these objects. In this chapter, we will focus on the functions, methods, and routines that are significant to matrices alone. We will disregard operations if their outputs have no translation into linear algebra equivalents. For a primer on `ndarray`, we recommend you to browse through Chapter 2, Top-Level SciPy of Learning SciPy for Numerical and Scientific Computing, Second Edition. For a quick review of Linear Algebra, we recommend Hoffman and Kunze, Linear Algebra 2nd Edition, Pearson, 1971.

### Constructing matrices in the ndarray class

We may create matrices from data as `ndarray` instances in three different ways: manually from standard input, by assigning to each entry a value from a function, or by retrieving the data from external files.

Constructor

Description

`numpy.array(object)`

Create a matrix from `object`

`numpy.diag(arr, k)`

Create diagonal matrix with entries of array `arr` on diagonal `k`

`numpy.fromfunction(function, shape)`

Create a matrix by executing a function over each coordinate

`numpy.fromfile(fname)`

Create a matrix from a text or binary file (basic)

`numpy.loadtxt(fname)`

Create a matrix from a text file (advanced)

Let us create some example matrices to illustrate some of the functions defined in the previous table. As before, we start an iPython session:

```In [1]: import numpy as np, matplotlib.pyplot as plt, \
...: scipy.linalg as spla, scipy.sparse as spsp, \
...: scipy.sparse.linalg as spspla
In [2]: A = np.array([[1,2],[4,16]]); \...: A
Out[2]:
array([[ 1,  2],
[ 4, 16]])
In [3]: B = np.fromfunction(lambda i,j: (i-1)*(j+1),
...:                     (3,2), dtype=int); \
...: print B
...:
[[-1 -2]
[ 0  0]
[ 1  2]]
In [4]: np.diag((1j,4))
Out[4]:
array([[ 0.+1.j,  0.+0.j],
[ 0.+0.j,  4.+0.j]])
```

Special matrices with predetermined zeros and ones can be constructed with the following functions:

Constructor

Description

`numpy.empty(shape)`

Array of a given shape, entries not initialized

`numpy.eye(N, M, k)`

2-D array with ones on the k-th diagonal, and zeros elsewhere

`numpy.identity(n)`

Identity array

`numpy.ones(shape)`

Array with all entries equal to one

`numpy.zeros(shape)`

Array with all entries equal to zero

`numpy.tri(N, M, k)`

Array with ones at and below the given diagonal, zeros otherwise

### Tip

All these constructions, except `numpy.tri`, have a companion function `xxx_like` that creates `ndarray` with the requested characteristics and with the same shape and data type as another source `ndarray` class:

```In [5]: np.empty_like(A)
Out[5]:
array([[140567774850560, 140567774850560],
[     4411734640, 562954363882576]])
```

Of notable importance are arrays constructed as numerical ranges.

Constructor

Description

`numpy.arange(stop)`

Evenly spaced values within an interval

`numpy.linspace(start, stop)`

Evenly spaced numbers over an interval

`numpy.logspace(start, stop)`

Evenly spaced numbers on a log scale

`numpy.meshgrid`

Coordinate matrices from two or more coordinate vectors

`numpy.mgrid`

`nd_grid` instance returning dense multi-dimensional `meshgrid`

`numpy.ogrid`

`nd_grid` instance returning open multi-dimensional `meshgrid`

Special matrices with numerous applications in linear algebra can be easily called from within NumPy and the module `scipy.linalg`.

Constructor

Description

`scipy.linalg.circulant(arr)`

Circulant matrix generated by 1-D array `arr`

`scipy.linalg.companion(arr)`

Companion matrix of polynomial with coefficients coded by `arr`

`scipy.linalg.hadamard(n)`

Sylvester's construction of a Hadamard matrix of size n × n. n must be a power of 2

`scipy.linalg.hankel(arr1, arr2)`

Hankel matrix with `arr1` as the first column and `arr2` as the last column

`scipy.linalg.hilbert(n)`

Hilbert matrix of size n × n

`scipy.linalg.invhilbert(n)`

The inverse of a Hilbert matrix of size n × n

`scipy.linalg.leslie(arr1, arr2)`

Leslie matrix with fecundity array `arr1` and survival coefficients `arr2`

`scipy.linalg.pascal(n)`

n × n truncations of the Pascal matrix of binomial coefficients

`scipy.linalg.toeplitz(arr1, arr2)`

Toeplitz array with first column `arr1` and first row `arr2`

`numpy.vander(arr)`

Van der Monde matrix of array `arr`

For instance, one fast way to obtain all binomial coefficients of orders up to a large number (the corresponding Pascal triangle) is by means of a precise Pascal matrix. The following example shows how to compute these coefficients up to order `13`:

```In [6]: print spla.pascal(13, kind='lower')
```

Besides these basic constructors, we can always stack arrays in different ways:

Constructor

Description

`numpy.concatenate((A1, A2, ...))`

Join matrices together

`numpy.hstack((A1, A2, ...))`

Stack matrices horizontally

`numpy.vstack((A1, A2, ...))`

Stack matrices vertically

`numpy.tile(A, reps)`

Repeat a matrix a certain number of times (given by `reps`)

`scipy.linalg.block_diag(A1,A2, ...)`

Create a block diagonal array

Let us observe some of these constructors in action:

```In [7]: np.tile(A, (2,3))   # 2 rows, 3 columns
Out[7]:
array([[ 1,  2,  1,  2,  1,  2],
[ 4, 16,  4, 16,  4, 16],
[ 1,  2,  1,  2,  1,  2],
[ 4, 16,  4, 16,  4, 16]])
In [8]: spla.block_diag(A,B)
Out[8]:
array([[ 1,  2,  0,  0],
[ 4, 16,  0,  0],
[ 0,  0, -1, -2],
[ 0,  0,  0,  0],
[ 0,  0,  1,  2]])
```

### Constructing matrices in the matrix class

For the `matrix` class, the usual way to create a matrix directly is to invoke either `numpy.mat` or `numpy.matrix`. Observe how much more comfortable is the syntax of `numpy.matrix` than that of `numpy.array`, in the creation of a matrix similar to `A`. With this syntax, different values separated by commas belong to the same row of the matrix. A semi-colon indicates a change of row. Notice the casting to the `matrix` class too!

```In [9]: C = np.matrix('1,2;4,16'); \
...: C
Out[9]:
matrix([[ 1,  2],
[ 4, 16]])
```

These two functions also transform any `ndarray` into `matrix`. There is a third function that accomplishes this task: `numpy.asmatrix`:

```In [10]: np.asmatrix(A)
Out[10]:
matrix([[ 1,  2],
[ 4, 16]])
```

For arrangements of matrices composed by blocks, besides the common stack operations for `ndarray` described before, we have the extremely convenient function `numpy.bmat`. Note the similarity with the syntax of `numpy.matrix`, particularly the use of commas to signify horizontal concatenation and semi-colons to signify vertical concatenation:

```In [11]: np.bmat('A;B')        In [12]: np.bmat('A,C;C,A')
Out[11]:                       Out[12]:
matrix([[ 1,  2],              matrix([[ 1,  2,  1,  2],
[ 4, 16],                      [ 4, 16,  4, 16],
[-1, -2],                      [ 1,  2,  1,  2],
[ 0,  0],                      [ 4, 16,  4, 16]])
[ 1,  2]])
```

### Constructing sparse matrices

There are seven different ways to input sparse matrices. Each format is designed to make a specific problem or operation more efficient. Let us go over them in detail:

Method

Name

Optimal use

`BSR`

Block Sparse Row

Efficient arithmetic, provided the matrix contains blocks.

`COO`

Coordinate

Fast and efficient construction format. Efficient methods to convert to the CSC and CSR formats.

`CSC`

Compressed Sparse Column

Efficient matrix arithmetic and column slicing. Relatively fast matrix-vector product.

`CSR`

Compressed Sparse Row

Efficient matrix arithmetic and row slicing. Fastest to perform matrix-vector products.

`DIA`

Diagonal storage

Efficient for construction and storage if the matrix contains long diagonals of non-zero entries.

`DOK`

Dictionary of keys

Efficient incremental construction and access of individual matrix entries.

`LIL`

Flexible slicing. Efficient for changes to matrix sparsity.

They can be populated in up to five ways, three of which are common to every sparse matrix format:

• They can cast to sparse any generic matrix. The `lil` format is the most effective with this method:

```In [13]: A_coo = spsp.coo_matrix(A); \
....: A_lil = spsp.lil_matrix(A)
```
• They can cast to a specific sparse format another sparse matrix in another sparse format:

```In [14]: A_csr = spsp.csr_matrix(A_coo)
```
• Empty sparse matrices of any shape can be constructed by indicating the shape and `dtype`:

```In [15]: M_bsr = spsp.bsr_matrix((100,100), dtype=int)
```

They all have several different extra input methods, each specific to their storage format.

• Fancy indexing: As we would do with any generic matrix. This is only possible with the LIL or DOK formats:

```In [16]: M_lil = spsp.lil_matrix((100,100), dtype=int)
In [17]: M_lil[25:75, 25:75] = 1
In [18]: M_bsr[25:75, 25:75] = 1
NotImplementedError    Traceback (most recent call last)
<ipython-input-18-d9fa1001cab8> in <module>()
----> 1 M_bsr[25:75, 25:75] = 1
[...]/scipy/sparse/bsr.pyc in __setitem__(self, key, val)
297
298     def __setitem__(self,key,val):
--> 299         raise NotImplementedError
300
301     ######################
NotImplementedError:
```
• Dictionary of keys: This input system is most effective when we create, update, or search each element one at a time. It is efficient only for the LIL and DOK formats:

```In [19]: M_dok = spsp.dok_matrix((100,100), dtype=int)
In [20]: position = lambda i, j: ((i<j) & ((i+j)%10==0))
In [21]: for i in range(100):
....:     for j in range(100):
....:         M_dok[i,j] = position(i,j)
....:
```
• Data, rows, and columns: This is common to four formats: BSR, COO, CSC, and CSR. This is the method of choice to import sparse matrices from the Matrix Market Exchange format, as illustrated at the beginning of the chapter.

### Tip

With the data, rows, and columns input method, it is a good idea to always include the option `shape` in the construction. In case this is not provided, the size of the matrix will be inferred from the largest coordinates from the rows and columns, resulting possibly in a matrix of a smaller size than required.

• Data, indices, and pointers: This is common to three formats: BSR, CSC, and CSR. It is the method of choice to import sparse matrices from the Rutherford-Boeing Exchange format.

### Note

The Rutherford-Boeing Exchange format is an updated version of the Harwell-Boeing format. It stores the matrix as three vectors: `pointers_v`, `indices_v`, and `data`. The row indices of the entries of the jth column are located in positions `pointers_v(j)` through `pointers_v(j+1)-1` of the vector `indices_v`. The corresponding values of the matrix are located at the same positions, in the vector data.

Let us show by example how to read an interesting matrix in the Rutherford-Boeing matrix exchange format, `Pajek/football`. This 35 × 35 matrix with 118 non-zero entries can be found in the collection at www.cise.ufl.edu/research/sparse/matrices/Pajek/football.html.

It is an adjacency matrix for a network of all the national football teams that attended the FIFA World Cup celebrated in France in 1998. Each node in the network represents one country (or national football team) and the links show which country exported players to another country.

This is a printout of the `football.rb` file:

The header of the file (the first four lines) contains important information:

• The first line provides us with the title of the matrix, `Pajek/football; 1998; L. Krempel; ed: V. Batagelj`, and a numerical key for identification purposes `MTRXID=1474`.

• The second line contains four integer values: `TOTCRD=12` (lines containing significant data after the header; see In [24]), `PTRCRD=2` (number of lines containing pointer data), `INDCRD=5` (number of lines containing indices data), and `VALCRD=2` (number of lines containing the non-zero values of the matrix). Note that it must be TOTCRD = PTRCRD + INDCRD + VALCRD.

• The third line indicates the matrix type `MXTYPE=(iua)`, which in this case stands for an integer matrix, unsymmetrical, compressed column form. It also indicates the number of rows and columns (`NROW=35`, `NCOL=35`), and the number of non-zero entries (`NNZERO=118`). The last entry is not used in the case of a compressed column form, and it is usually set to zero.

• The fourth column contains the Fortran formats for the data in the following columns. `PTRFMT=(20I4)` for the pointers, `INDFMT=(26I3)` for the indices, and `VALFMT=(26I3)` for the non-zero values.

We proceed to opening the file for reading, storing each line after the header in a Python list, and extracting from the relevant lines of the file, the data we require to populate the vectors `indptr`, `indices`, and `data`. We finish by creating the corresponding sparse matrix called `football` in the CSR format, with the `data`, `indices`, `pointers` method:

```In [22]: f = open("football.rb", 'r'); \
....: football_list = list(f); \
....: f.close()
In [23]: football_data = np.array([])
In [24]: for line in range(4, 4+12):
....:     newdata = np.fromstring(football_list[line], sep=" ")
....:     football_data = np.append(football_data, newdata)
....:
In [25]: indptr = football_data[:35+1] - 1; \
....: indices = football_data[35+1:35+1+118] - 1; \
....: data = football_data[35+1+118:]
In [26]: football = spsp.csr_matrix((data, indices, indptr),
....:                            shape=(35,35))
```

At this point, it is possible to visualize the network with its associated graph, with the help of a Python module called `networkx`. We obtain the following diagram depicting as nodes the different countries. Each arrow between the nodes indicates the fact that the originating country has exported players to the receiving country:

### Note

`networkx` is a Python module to deal with complex networks. For more information, visit their Github project pages at networkx.github.io.

One way to accomplish this task is as follows:

```In [27]: import networkx
In [28]: G = networkx.DiGraph(football)
In [29]: f = open("football_nodename.txt"); \
....: m = list(f); \
....: f.close()
In [30]: def rename(x): return m[x]
In [31]: G = networkx.relabel_nodes(G, rename)
In [32]: pos = networkx.spring_layout(G)
In [33]: networkx.draw_networkx(G, pos, alpha=0.2, node_color='w',
....:                        edge_color='b')
```

The module `scipy.sparse` borrows from NumPy some interesting concepts to create constructors and special matrices:

Constructor

Description

`scipy.sparse.diags(diagonals, offsets)`

Sparse matrix from diagonals

`scipy.sparse.rand(m, n, density)`

Random sparse matrix of prescribed density

`scipy.sparse.eye(m)`

Sparse matrix with ones in the main diagonal

`scipy.sparse.identity(n)`

Identity sparse matrix of size n × n

Both functions `diags` and `rand` deserve examples to show their syntax. We will start with a sparse matrix of size 14 × 14 with two diagonals: the main diagonal contains 1s, and the diagonal below contains 2s. We also create a random matrix with the function `scipy.sparse.rand`. This matrix has size 5 × 5, with 25 percent non-zero elements (`density=0.25`), and is crafted in the LIL format:

```In [34]: diagonals = [[1]*14, [2]*13]
In [35]: print spsp.diags(diagonals, [0,-1]).todense()
[[ 1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
[ 2.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
[ 0.  2.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
[ 0.  0.  2.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
[ 0.  0.  0.  2.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
[ 0.  0.  0.  0.  2.  1.  0.  0.  0.  0.  0.  0.  0.  0.]
[ 0.  0.  0.  0.  0.  2.  1.  0.  0.  0.  0.  0.  0.  0.]
[ 0.  0.  0.  0.  0.  0.  2.  1.  0.  0.  0.  0.  0.  0.]
[ 0.  0.  0.  0.  0.  0.  0.  2.  1.  0.  0.  0.  0.  0.]
[ 0.  0.  0.  0.  0.  0.  0.  0.  2.  1.  0.  0.  0.  0.]
[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  2.  1.  0.  0.  0.]
[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  2.  1.  0.  0.]
[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  2.  1.  0.]
[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  2.  1.]]
In [36]: S_25_lil = spsp.rand(5, 5, density=0.25, format='lil')
In [37]: S_25_lil
Out[37]:
<5x5 sparse matrix of type '<type 'numpy.float64'>'
with 6 stored elements in LInked List format>
In [38]: print S_25_lil
(0, 0)    0.186663044982
(1, 0)    0.127636181284
(1, 4)    0.918284870518
(3, 2)    0.458768884701
(3, 3)    0.533573291684
(4, 3)    0.908751420065
In [39]: print S_25_lil.todense()
[[ 0.18666304  0.          0.          0.          0.        ]
[ 0.12763618  0.          0.          0.          0.91828487]
[ 0.          0.          0.          0.          0.        ]
[ 0.          0.          0.45876888  0.53357329  0.        ]
[ 0.          0.          0.          0.90875142  0.        ]]
```

Similar to the way we combined `ndarray` instances, we have some clever ways to combine sparse matrices to construct more complex objects:

Constructor

Description

`scipy.sparse.bmat(blocks)`

Sparse matrix from sparse sub-blocks

`scipy.sparse.hstack(blocks)`

Stack sparse matrices horizontally

`scipy.sparse.vstack(blocks)`

Stack sparse matrices vertically

### Linear operators

A linear operator is basically a function that takes as input a column vector and outputs another column vector, by left multiplication of the input with a matrix. Although technically, we could represent these objects just by handling the corresponding matrix, there are better ways to do this.

Constructor

Description

`scipy.sparse.linalg.LinearOperator(shape, matvec)`

Common interface for performing matrix vector products

`scipy.sparse.linalg.aslinearoperator(A)`

Return `A` as `LinearOperator`

In the `scipy.sparse.linalg` module, we have a common interface that handles these objects: the `LinearOperator` class. This class has only the following two attributes and three methods:

• `shape`: The shape of the representing matrix

• `dtype`: The data type of the matrix

• `matvec`: To perform multiplication of a matrix with a vector

• `rmatvec`: To perform multiplication by the conjugate transpose of a matrix with a vector

• `matmat`: To perform multiplication of a matrix with another matrix

Its usage is best explained through an example. Consider two functions that take vectors of size 3, and output vectors of size 4, by left multiplication with two respective matrices of size 4 × 3. We could very well define these functions with lambda predicates:

```In [40]: H1 = np.matrix("1,3,5; 2,4,6; 6,4,2; 5,3,1"); \
....: H2 = np.matrix("1,2,3; 1,3,2; 2,1,3; 2,3,1")
In [41]: L1 = lambda x: H1.dot(x); \
....: L2 = lambda x: H2.dot(x)
In [42]: print L1(np.ones(3))
[[  9.  12.  12.   9.]]
In [43]: print L2(np.tri(3,3))
[[ 6.  5.  3.]
[ 6.  5.  2.]
[ 6.  4.  3.]
[ 6.  4.  1.]]
```

Now, one issue arises when we try to add/subtract these two functions, or multiply any of them by a scalar. Technically, it should be as easy as adding/subtracting the corresponding matrices, or multiplying them by any number, and then performing the required left multiplication again. But that is not the case.

For instance, we would like to write `(L1+L2)(v)` instead of `L1(v) + L2(v)`. Unfortunately, doing so will raise an error:

```TypeError: unsupported operand type(s) for +: 'function' and
'function'
```

Instead, we may instantiate the corresponding linear operators and manipulate them at will, as follows:

```In [44]: Lo1 = spspla.aslinearoperator(H1); \
....: Lo2 = spspla.aslinearoperator(H2)
In [45]: Lo1 - 6 * Lo2
Out[45]: <4x3 _SumLinearOperator with dtype=float64>
In [46]: print Lo1 * np.ones(3)
[  9.  12.  12.  9.]
In [47]: print (Lo1-6*Lo2) * np.tri(3,3)
[[-27. -22. -13.]
[-24. -20.  -6.]
[-24. -18. -16.]
[-27. -20.  -5.]]
```

Linear operators are a great advantage when the amount of information needed to describe the product with the related matrix is less than the amount of memory needed to store the non-zero elements of the matrix.

For instance, a permutation matrix is a square binary matrix (ones and zeros) that has exactly one entry in each row and each column. Consider a large permutation matrix, say 1024 × 1024, formed by four blocks of size 512 × 512: a zero block followed horizontally by an identity block, on top of an identity block followed horizontally by another zero block. We may store this matrix in three different ways:

```In [47]: P_sparse = spsp.diags([[1]*512, [1]*512], [512,-512], \
....:                       dtype=int)
In [48]: P_dense = P_sparse.todense()
In [49]: mv = lambda v: np.roll(v, len(v)/2)
In [50]: P_lo = spspla.LinearOperator((1024,1024), matvec=mv, \
....:                              matmat=mv, dtype=int)
```

In the sparse case, `P_sparse`, we may think of this as the storage of just 1024 integer numbers. In the dense case, `P_dense`, we are technically storing 1048576 integer values. In the case of the linear operator, it actually looks like we are not storing anything! The function `mv` that indicates how to perform the multiplications has a much smaller footprint than any of the related matrices. This is also reflected in the time of execution of the multiplications with these objects:

```In [51]: %timeit P_sparse * np.ones(1024)
10000 loops, best of 3: 29.7 µs per loop
In [52]: %timeit P_dense.dot(np.ones(1024))
100 loops, best of 3: 6.07 ms per loop
In [53]: %timeit P_lo * np.ones(1024)
10000 loops, best of 3: 25.4 µs per loop
```