Programming with Python

Stefan Güttel, guettel.com

Lab classes: Modules

Problem 1. Write a module symtoep for creating various symmetric Toeplitz matrices with the following functions:

  • general(L) that returns a symmetric Toeplitz matrix with the first row equal to the elements of the list L.

  • tridiagonal(n, d, sd) that returns a tridiagonal symmetric Toeplitz matrix of order n, with the float number d on the diagonal and the float number sd right below and above it (the rest of the elements are zero).

  • filled(n, d, nd) that returns a symmetric Toeplitz matrix of order n, with the float number d on the diagonal and all the non-diagonal elements equal to nd.

  • menu() that displays a choice

     Load a symmetric Toeplitz matrix:
   1. General symmetric Toeplitz matrix
   2. Tridiagonal symmetric Toeplitz matrix
   3. Filled symmetric Toeplitz matrix
 Your choice (1-3):

    and asks the user to choose. Then,

    • if the user chooses "1", the function asks for a list L of numbers (preferably as a string of comma-separated floats, but you can use some other method as well), and then returns general(L),
    • if the user chooses "2", the function asks for an integer n and floats d and sd, and returns tridiagonal(n, d, sd).
    • if the user chooses "3", the function asks for an integer n and floats d and nd, and returns filled(n, d, sd).
    • if the user chooses anything else, the choice is redisplayed and user is asked to choose again.

When run like a program, this module should do nothing.

Toeplitz matrices are those that have constant elements on all diagonals.

  • An example of a symmetric Toeplitz matrix of order $5$: $$\begin{bmatrix} 1 & 2 & 3 & 4 & 5 \\ 2 & 1 & 2 & 3 & 4 \\ 3 & 2 & 1 & 2 & 3 \\ 4 & 3 & 2 & 1 & 2 \\ 5 & 4 & 3 & 2 & 1 \end{bmatrix}$$
  • An example of a tridiagonal (all elements other than those on the main and its neighbour diagonals are zero) symmetric Toeplitz matrix of order $5$: $$\begin{bmatrix} 1 & 2 & 0 & 0 & 0 \\ 2 & 1 & 2 & 0 & 0 \\ 0 & 2 & 1 & 2 & 0 \\ 0 & 0 & 2 & 1 & 2 \\ 0 & 0 & 0 & 2 & 1 \end{bmatrix}$$
  • An example of a "filled" symmetric Toeplitz matrix of order $5$: $$\begin{bmatrix} 1 & 2 & 2 & 2 & 2 \\ 2 & 1 & 2 & 2 & 2 \\ 2 & 2 & 1 & 2 & 2 \\ 2 & 2 & 2 & 1 & 2 \\ 2 & 2 & 2 & 2 & 1 \end{bmatrix}$$

All the functions have to return the matrices as NumPy's matrix type, which is easily created from a list or any other two-dimensional array-like type, using the function numpy.matrix. For example, the first of the three matrices shown above can be created as follows:

In [1]:
import numpy as np
L = list(range(1,6))
res = list()
for i in range(5):
    row = list()
    for j in range(5):
        row.append(L[abs(i-j)])
    res.append(row)
res = np.matrix(res)
print(res)

[[1 2 3 4 5]
 [2 1 2 3 4]
 [3 2 1 2 3]
 [4 3 2 1 2]
 [5 4 3 2 1]]

It can be done more easily using a generator expression:

In [2]:
import numpy as np
L = list(range(1,6))
res = np.matrix([[L[abs(i-j)] for i in range(5)] for j in range(5)])
print(res)

[[1 2 3 4 5]
 [2 1 2 3 4]
 [3 2 1 2 3]
 [4 3 2 1 2]
 [5 4 3 2 1]]

You may implement the functions tridiagonal(n, d, sd) and filled(n, d, nd) either on their own or simply by creating L and calling general(L). The former approach is prefered, in which case numpy.zeros, numpy.ones, numpy.fill_diagonal, and numpy.tri may help you, along with the matrix addition (done with the usual plus + operator) and multiplication with a constant (done with the usual multiplication * operator). These functions are much faster than setting the elements index by index, but you are free to work without them.

Problem 2. Using the module symtoep, write a program that loads a matrix (one of the 3 supported types), and then prints it and its eigenvalues.

To compute the eigenvalues of a symmetric matrix, use scipy.linalg.eigvalsh. If, for some reason, SciPy doesn't work, use numpy.linalg.eigvalsh.

Ideally, make your program run properly with either (but this is not a must). For example,

try:
    import scipy.linalg as la
except ImportError:
    import numpy.linalg as la
...
print(la.eigvalsh(A))

Problem 3. Using the module symtoep, write a program that loads a matrix (one of the 3 supported types), and then prints it and a message explaining if it is positive semidefinite or not.

To test positive semidefiniteness, try to compute the Cholesky factorization, using either scipy.linalg.cholesky or numpy.linalg.cholesky. If the computation fails (i.e., a numpy.linalg.linalg.LinAlgError exception is raised), the matrix is not positive semidefinite. Otherwise, it is.

The same remarks from Problem 2 regarding the import of SciPy/NumPy modules apply.