Python:Linear Algebra

General

Solving

To solve Ax=b using linear algebra, be sure that A is a 2D array. b can either be 1D or 2D -- and in fact if 2D it can be a row or a column! Some math packages that solve linear algebra problems would require that b be a 2D column, but not Python. The result x will be the same shape and size as b (that is, 1D, 2D row, or 2D column). Here is an example with b as a 2D column:

A = np.array([[1, 1], [1, -1]])
b = np.array([[3], [4]])
soln = np.linalg.solve(A, b)

for example will create a variable called soln that is:

array([[ 3.5],
       [-0.5]])

Printing solution

Python's format command is picky and will not take arrays. To print out the solutions to the above, for instance, you would need:

print('x: {:f}, y: {:f}'.format(soln[0][0], soln[1][0]))

or perhaps more clearly:

soln_vec = soln[:,0]
print('x: {:f}, y: {:f}'.format(soln_vec[0], soln_vec[1]))

Note that soln[0] will give you an array containing the $$x$$ value, not just the $$x$$ value, and similarly soln[1] will give you an array containing the $$y$$ value, not just the $$y$$ value, so while:

print('x: {}, y: {}'.format(soln[0], soln[1]))

works, it is printing out arrays:

x: [-4.], y: [4.5]

That means you cannot try to format the arrays as if they were floats:

print('x: {:f}, y: {:f}'.format(soln[0], soln[1]))

will give an error:

TypeError: unsupported format string passed to numpy.ndarray.__format__

Norms

Vector norms and matrix norms are different ways of quantifying the "size" of an array, not in terms of the number of rows or columns but in terms of the values of the entries themselves.

1D Arrays

For 1D arrays only (i.e. a column vector or row vector), the $$p$$ norm for a vector $$x$$, $$||x||_p$$, is defined as:

\( ||x||_p=\left(\sum_k |x_k|^p\right)^{1/p} \)

Typical values of $$p$$ include 1, 2, and $$\infty$$; for the latter, the $$p$$ norm is defined as the largest absolute value in the array (regardless of whether it is repeated). The 2 norm is also known as the Euclidean norm and thus may also be denoted $$||x||_e$$

2D Arrays

For 2D arrays, there are four common norms. Assuming some two-dimensional array $$A$$, they are:

• Matrix 1 norm $$||A||_1$$: The largest 1 norm of the columns of $$A$$
• Matrix $$\infty$$ norm $$||A||_{\infty}$$: The largest 1 norm of the rows of $$A$$ (note that the matrix $$\infty$$ norm is defined based on the vector 1 norm
• Matrix Frobenius norm $$||A||_{f}$$: The square root of the sum of the squares of the absolute values of the entries in the matrix; \(||A||_f=\sqrt{\sum_i\sum_j|a_{ij}|^2}\)where $$a_{i,j}$$ is the entry on row $$i$$ and column $$j$$ of matrix $$A$$. The absolute value is needed in case entries are complex numbers. This is basically the Euclidean norm for a matrix.
• Matrix spectral norm or 2 norm $$||A||_2$$: the square root of the largest eigenvalue of $$A^{T}A$$. Note that this is calculated in a vastly different way from the 2 norm of a 1D array!

Condition numbers

As noted in Chapra 11.2.2, the base-10 logarithm gives an estimate for how many digits of precision are lost between the number of digits in the coefficients and the number of digits in the solution. Condition numbers based on the 2-norm may be calculated in Python using:

np.linalg.cond(A) # default is based on 2-norm

For information on using other norms to calculate condition numbers, see

help(np.linalg.cond)

and specifically information about the kwarg p. Note that you must use np.inf, not just inf, for the infinity norm.

Sweeping a Parameter

If you have a system where the coefficients change as a function of some parameter, you will generally need to use a loop to solve and store the solutions. If you have a system where the forcing function (right-side vector) changes, you may be able to solve all at once but generally a loop is the way to go. The following shows example code for sweeping through a parameter, storing values, and then plotting them:

Changing coefficient matrix

Equations

For this example, the equations are:

\( \begin{align} mx-y&=4\\ x+y&=3 \end{align} \)

which means a matrix-based representation is:

\( \begin{align} \begin{bmatrix} m & -1\\1 & 1 \end{bmatrix} \begin{bmatrix} x\\y \end{bmatrix} &= \begin{bmatrix}4\\3 \end{bmatrix} \end{align} \)

The determinant for the coefficient matrix of this system is \(m+1\) meaning there should be a unique solution for all values of \(m\) other than -1. The code is going to sweep through 50 values of \(m\) between 0 and 15.

Alternate Loop Structures

There are at least three ways to control the loop in the code above and at least two different ways to store the values. As shown, the loop uses a variable k and assigns it to the value in range(npts), which will contain integers from 0 to npts-1 by 1.

Loop Over Slope Values With External Counter

The loop could be made to iterate through the different values of marray, though that would require an external counter:

counter = 0
for slope in marray:
    A = np.array([[slope, -1], [1, 1]])
    b = np.array([[4], [3]])
    soln = np.linalg.solve(A, b)
    x[counter] = soln[0][0]
    y[counter] = soln[1][0]
    counter += 1

In this case, slope will iterate through each of the slope values in marray. The slope in the A line could be replaced with marray[counter], but that would defeat the purpose of iterating through m

Loop Over Enumerated Slope Values

In a "best of both worlds" situation, the loop could be made to iterate through an enumeration of the different values of m. This way, the loop variables would keep track of both the loop count and the specific slope:

for k, slope in enumerate(m):
    A = np.array([[slope, -1], [1, 1]])
    b = np.array([[4], [3]])
    soln = np.linalg.solve(A, b)
    x[k] = soln[0][0]
    y[k] = soln[1][0]

In this case, the slope in the A line could be replaced with m[k], but that would defeat the purpose of iterating through the values of m!

Using Lists

In each of the above cases, the x and y values could be stored in lists; new values would be appended to the list and that does not require an external counter. The initialization would also be different:

npts = 50
m = np.linspace(0, 5, npts)
xL = []
yL = []

for slope in m:
    A = np.array([[slope, -1], [1, 1]])
    b = np.array([[4], [3]])
    soln = np.linalg.solve(A, b)
    xL.append(soln[0][0])
    yL.append(soln[1][0])

x = np.array(xL)
y = np.array(yL)

Changing solution vector

Equations

For this example, the equations are:

\( \begin{align} x-y&=p\\ x+y&=3 \end{align} \)

which means a matrix-based representation is:

\( \begin{align} \begin{bmatrix} 1 & -1\\1 & 1 \end{bmatrix} \begin{bmatrix} x\\y \end{bmatrix} &= \begin{bmatrix}p\\3 \end{bmatrix} \end{align} \)

The determinant for the coefficient matrix of this system is 2 meaning there should always be a unique solution. The code is going to sweep through 75 values of \(p\) between -5 and 10.

Note: the condition number for this system will remain constant since the $$A$$ matrix does not change.

Multiple solution vectors simultaneously

This method is not recommended for people with limited experience with linear algebra.

Equations

For this example, the equations are:

\( \begin{align} x-y&=p\\ x+y&=3 \end{align} \)

which means a matrix-based representation is:

\( \begin{align} \begin{bmatrix} 1 & -1\\1 & 1 \end{bmatrix} \begin{bmatrix} x\\y \end{bmatrix} &= \begin{bmatrix}p\\3 \end{bmatrix} \end{align} \)

The determinant for the coefficient matrix of this system is 2 meaning there should always be a unique solution. The code is going to solve the system for 75 values of \(p\) between -5 and 10 by setting up a 75-column matrix of solution vectors and then extracting the first row of solutions for \(xa\) and the second row for \(ya\). Note that unlike the above two examples where \(x\) and \(y\) were lists, \(xa\) and \(ya\) are arrays.

Code

import numpy as np
import matplotlib.pyplot as plt

p = np.linspace(-5, 10, 75)
rhs = np.block([[p], [3 + 0 * p]]) # note use of 0*p to get array of correct size!
all_soln = np.linalg.solve(A, rhs)
xa = all_soln[0][:]
ya = all_soln[1][:]

Sweeping Two Parameters

If there are two quantities that change in your system, you can generate a surface of solutions by using a double for-loop to solve your system and store individual solutions in an appropriately-sized array, then use 3D graphics to produce a surface.