Difference between revisions of "Python:Linear Algebra"

From PrattWiki
Jump to navigation Jump to search
(Sweeping Two Parameters)
 
(27 intermediate revisions by the same user not shown)
Line 1: Line 1:
 +
== 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:
 +
<source lang=python>
 +
A = np.array([[1, 1], [1, -1]])
 +
b = np.array([[3], [4]])
 +
soln = np.linalg.solve(A, b)
 +
</source>
 +
for example will create a variable called <code>soln</code> that is:
 +
<source lang=python>
 +
array([[ 3.5],
 +
      [-0.5]])
 +
</source>
 +
 +
===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:
 +
<syntaxhighlight lang=python>
 +
print('x: {:f}, y: {:f}'.format(soln[0][0], soln[1][0]))
 +
</syntaxhighlight >
 +
or perhaps more clearly:
 +
<syntaxhighlight lang=python>
 +
soln_vec = soln[:,0]
 +
print('x: {:f}, y: {:f}'.format(soln_vec[0], soln_vec[1]))
 +
</syntaxhighlight >
 +
Note that <code>soln[0]</code> will give you an ''array'' containing the $$x$$ value, not just the $$x$$ value, and similarly <code>soln[1]</code> will give you an ''array'' containing the $$y$$ value, not just the $$y$$ value, so while:
 +
<syntaxhighlight lang=python>
 +
print('x: {}, y: {}'.format(soln[0], soln[1]))
 +
</syntaxhighlight>
 +
works, it is printing out arrays:
 +
<source>
 +
x: [-4.], y: [4.5]
 +
</source>
 +
 +
That means you cannot try to format the arrays as if they were floats:
 +
<syntaxhighlight  lang=python>
 +
print('x: {:f}, y: {:f}'.format(soln[0], soln[1]))
 +
</syntaxhighlight >
 +
will give an error:
 +
<syntaxhighlight >
 +
TypeError: unsupported format string passed to numpy.ndarray.__format__
 +
</syntaxhighlight >
 +
 +
=== 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:
 +
<center><math>
 +
||x||_p=\left(\sum_k |x_k|^p\right)^{1/p}
 +
</math></center>
 +
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; <center><math>||A||_f=\sqrt{\sum_i\sum_j|a_{ij}|^2}</math></center>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:
 +
<source lang=python>
 +
np.linalg.cond(A) # default is based on 2-norm
 +
</source>
 +
For information on using other norms to calculate condition numbers, see
 +
<source lang=python>
 +
help(np.linalg.cond)
 +
</source>
 +
and specifically information about the kwarg <code>p</code>.  Note that you must use <code>np.inf</code>, not just inf, for the infinity norm.
 +
 
== Sweeping a Parameter ==
 
== 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:
 
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:
Line 21: Line 92:
 
\end{align}
 
\end{align}
 
</math></center>
 
</math></center>
The determinant for the coefficient matrix of this system is <math>m+1</math> meaning there should be a unique solution for all values of <math>m</math> other than -1.  The code is going to sweep through 50 values of <math>m</math> between 0 and 5.
+
The determinant for the coefficient matrix of this system is <math>m+1</math> meaning there should be a unique solution for all values of <math>m</math> other than -1.  The code is going to sweep through 50 values of <math>m</math> between 0 and 15.
  
 
==== Code ====
 
==== Code ====
 +
<html>
 +
<iframe src="https://trinket.io/embed/python3/b45f48df59?start=result" width="100%" height="600" frameborder="0" marginwidth="0" marginheight="0" allowfullscreen></iframe>
 +
</html>
 +
<!--
 
[[File:Lin_alg_demo.png|thumb|Graph of solutions with parameter sweep]]
 
[[File:Lin_alg_demo.png|thumb|Graph of solutions with parameter sweep]]
 
<source lang=python>
 
<source lang=python>
Line 47: Line 122:
 
plt.legend()
 
plt.legend()
 
</source>
 
</source>
 +
-->
 +
 +
==== 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 <code>k</code> and assigns it to the value in <code>range(npts)</code>, which will contain integers from 0 to <code>npts-1</code> by 1. 
 +
===== Loop Over Slope Values With External Counter =====
 +
The loop could be made to iterate through the different values of <code>marray</code>, though that would require an external counter:
 +
 +
<syntaxhighlight lang=python>
 +
counter = 0
 +
for m in marray:
 +
    A = np.array([[m, -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
 +
</syntaxhighlight>  In this case, <code>m</code> will iterate through each of the slope values in <code>marray</code>.  The <code>m</code> in the <code>A</code> line could be replaced with <code>marray[counter]</code>, but that would defeat the purpose of iterating through <code>marray</code>
 +
 +
===== 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 <code>marray</code>.  This way, the loop variables would keep track of both the loop count and the specific slope:
 +
<syntaxhighlight lang=python>
 +
for k, m in enumerate(marray):
 +
    A = np.array([[m, -1], [1, 1]])
 +
    b = np.array([[4], [3]])
 +
    soln = np.linalg.solve(A, b)
 +
    x[k] = soln[0][0]
 +
    y[k] = soln[1][0]
 +
</syntaxhighlight>  In this case, the <code>m</code> in the <code>A</code> line could be replaced with <code>marray[k]</code>, but that would defeat the purpose of iterating through the values of <code>marray</code>!
 +
 +
===== Using Lists =====
 +
In each of the above cases, the <code>x</code> and <code>y</code> 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:
 +
<syntaxhighlight lang=python>
 +
npts = 50
 +
marray = np.linspace(0, 5, npts)
 +
xL = []
 +
yL = []
 +
 +
for m in marray:
 +
    A = np.array([[m, -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)
 +
</syntaxhighlight >
  
 
=== Changing solution vector ===
 
=== Changing solution vector ===
Line 69: Line 191:
 
</math></center>
 
</math></center>
 
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 <math>p</math> between -5 and 10.
 
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 <math>p</math> between -5 and 10.
 +
 +
Note: the condition number for this system will remain constant since the $$A$$ matrix does not change.
  
 
==== Code ====
 
==== Code ====
 +
<html>
 +
<iframe src="https://trinket.io/embed/python3/13a111c7e8?start=result" width="100%" height="600" frameborder="0" marginwidth="0" marginheight="0" allowfullscreen></iframe>
 +
</html>
 +
<!--
 
[[File:Lin_alg_demo_2.png|thumb|Graph of solutions with parameter sweep]]
 
[[File:Lin_alg_demo_2.png|thumb|Graph of solutions with parameter sweep]]
 
<source lang=python>
 
<source lang=python>
Line 94: Line 222:
 
plt.legend()
 
plt.legend()
 
</source>
 
</source>
 +
-->
  
 
=== Multiple solution vectors simultaneously ===
 
=== Multiple solution vectors simultaneously ===
Line 116: Line 245:
 
\end{align}
 
\end{align}
 
</math></center>
 
</math></center>
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 <math>p</math> between -5 and 10 by setting up a 75-column matrix of solution vectors and then extracting the first row of solutions for <math>xa</math> and the second row for <math>ya</math>. Note that unlike the above two examples where <math>x</math> and <math>x</math> were lists, <math>xa</math> and <math>ya</math> are arrays.
+
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 <math>p</math> between -5 and 10 by setting up a 75-column matrix of solution vectors and then extracting the first row of solutions for <math>xa</math> and the second row for <math>ya</math>. Note that unlike the above two examples where <math>x</math> and <math>y</math> were lists, <math>xa</math> and <math>ya</math> are arrays.
  
 
==== Code ====
 
==== Code ====
Line 123: Line 252:
 
import matplotlib.pyplot as plt
 
import matplotlib.pyplot as plt
  
p = np.linspace(-5, 10, 75)
+
parray = np.linspace(-5, 10, 75)
rhs = np.block([[p], [3 + 0 * p]]) # note use of 0*p to get array of correct size!
+
rhs = np.block([[parray], [3 + 0 * parray]]) # note use of 0*parray to get array of correct size!
 
all_soln = np.linalg.solve(A, rhs)
 
all_soln = np.linalg.solve(A, rhs)
 
xa = all_soln[0][:]
 
xa = all_soln[0][:]
 
ya = all_soln[1][:]
 
ya = all_soln[1][:]
 
</source>
 
</source>
 +
 +
== 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. 
 +
 +
For example, for a relatively simple electric circuit made up of a voltage source $$v_S$$ and four resistors, Kirchhoff's Current Law can give us a system of two equations to solve for two unknown node voltages:
 +
 +
$$
 +
\begin{align*}
 +
\frac{v_x-v_s}{R_1}+\frac{v_x}{R_2}+\frac{v_x-v_y}{R_3}=0\\
 +
\frac{v_y-v_x}{R_3}+\frac{v_y}{R_4}=0
 +
\end{align*}
 +
$$
 +
 +
Written as a matrix with $$v_x$$ and $$v_y$$ as the unknowns, this becomes:
 +
 +
$$
 +
\begin{align*}
 +
\begin{bmatrix}
 +
\frac{1}{R_1}+\frac{1}{R_2}+\frac{1}{R_3} & -\frac{1}{R_3}\\
 +
-\frac{1}{R_3} & \frac{1}{R_3}+\frac{1}{R_4}
 +
\end{bmatrix}
 +
\begin{bmatrix}
 +
v_x\\v_y
 +
\end{bmatrix}&=
 +
\begin{bmatrix}
 +
\frac{v_s}{R_1}\\0\end{bmatrix}
 +
\end{align*}
 +
$$
 +
 +
To look at how $$v_x$$ and $$v_y$$ change as $$R_1$$ and $$R_4$$ change, you can create two matrices of $$R_1$$ and $$R_4$$ values and then use a double for loop to traverse the matrices:
 +
<HTML>
 +
<iframe src="https://trinket.io/embed/python3/df7b55f65c" width="100%" height="600" frameborder="0" marginwidth="0" marginheight="0" allowfullscreen></iframe>
 +
</html>

Latest revision as of 20:24, 3 March 2023

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.

Code

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 m in marray:
    A = np.array([[m, -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, m will iterate through each of the slope values in marray. The m in the A line could be replaced with marray[counter], but that would defeat the purpose of iterating through marray

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 marray. This way, the loop variables would keep track of both the loop count and the specific slope:

for k, m in enumerate(marray):
    A = np.array([[m, -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 m in the A line could be replaced with marray[k], but that would defeat the purpose of iterating through the values of marray!

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
marray = np.linspace(0, 5, npts)
xL = []
yL = []

for m in marray:
    A = np.array([[m, -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.

Code

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

parray = np.linspace(-5, 10, 75)
rhs = np.block([[parray], [3 + 0 * parray]]) # note use of 0*parray 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.

For example, for a relatively simple electric circuit made up of a voltage source $$v_S$$ and four resistors, Kirchhoff's Current Law can give us a system of two equations to solve for two unknown node voltages:

$$ \begin{align*} \frac{v_x-v_s}{R_1}+\frac{v_x}{R_2}+\frac{v_x-v_y}{R_3}=0\\ \frac{v_y-v_x}{R_3}+\frac{v_y}{R_4}=0 \end{align*} $$

Written as a matrix with $$v_x$$ and $$v_y$$ as the unknowns, this becomes:

$$ \begin{align*} \begin{bmatrix} \frac{1}{R_1}+\frac{1}{R_2}+\frac{1}{R_3} & -\frac{1}{R_3}\\ -\frac{1}{R_3} & \frac{1}{R_3}+\frac{1}{R_4} \end{bmatrix} \begin{bmatrix} v_x\\v_y \end{bmatrix}&= \begin{bmatrix} \frac{v_s}{R_1}\\0\end{bmatrix} \end{align*} $$

To look at how $$v_x$$ and $$v_y$$ change as $$R_1$$ and $$R_4$$ change, you can create two matrices of $$R_1$$ and $$R_4$$ values and then use a double for loop to traverse the matrices: