13. Linear Algebra

In [1]:

import math
import numpy as np
%pylab inline

%pylab is deprecated, use %matplotlib inline and import the required libraries.
Populating the interactive namespace from numpy and matplotlib

Linear Algebra¶

General Manipulations of Matrices¶

You can use regular NumPy arrays or you can use a special matrix class that offers some short cuts.

In [2]:

a = np.array([[1.0, 2.0], [3.0, 4.0]])

In [3]:

print(a)

[[1. 2.]
 [3. 4.]]

In [4]:

print(a.transpose())
print(a.T)

[[1. 3.]
 [2. 4.]]
[[1. 3.]
 [2. 4.]]

In [5]:

ainv = np.linalg.inv(a)
print(ainv)

[[-2.   1. ]
 [ 1.5 -0.5]]

In [6]:

print(np.dot(a, ainv))

[[1.0000000e+00 0.0000000e+00]
 [8.8817842e-16 1.0000000e+00]]

In [7]:

a @ ainv

Out[7]:

array([[1.0000000e+00, 0.0000000e+00],
       [8.8817842e-16, 1.0000000e+00]])

The eye() function will generate an identity matrix (as will the identity()).

In [8]:

print(np.eye(2))
print(np.identity(2))

[[1. 0.]
 [0. 1.]]
[[1. 0.]
 [0. 1.]]

We can solve Ax = b.

In [9]:

b = np.array([5, 7])
x = np.linalg.solve(a, b)
print(x)

[-3.  4.]

The matrix Class¶

In [10]:

A = np.matrix('1.0 2.0; 3.0 4.0')
print(A)
print(A.T)

[[1. 2.]
 [3. 4.]]
[[1. 3.]
 [2. 4.]]

In [11]:

X = np.matrix('5.0 7.0')
Y = X.T

print(A*Y)

[[19.]
 [43.]]

In [12]:

print(A.I*Y)

[[-3.]
 [ 4.]]

Tridiagonal Matrix Solve¶

Here we'll solve the equation:

$$f^{\prime\prime} = g(x)$$

with $g(x) = sin(x)$, and the domain $x \in [0, 2\pi]$, with boundary conditions $f(0) = f(2\pi) = 0$. The solution is simply $f(x) = -sin(x)$.

We'll use a grid of $N$ points with $x_0$ on the left boundary and $x_{N-1}$ on the right boundary.

We difference our equation as:

$$f_{i+1} - 2 f_i + f_{i-1} = \Delta x^2 g_i$$

We keep the boundary points fixed, so we only need to solve for the $N-2$ interior points. Near the boundaries, our difference is: $$f_2 - 2 f_1 = \Delta x^2 g_1$$

and

$$-2f_{N-1} + f_{N-2} = \Delta x^2 g_{N-1}$$

.

We can write the system of equations for solving for the $N-2$ interior points as:

\begin{equation} A = \left ( \begin{array}{ccccccc} -2 & 1 & & & & & \newline 1 & -2 & 1 & & & & \newline & 1 & -2 & 1 & & & \newline & & \ddots & \ddots & \ddots & & \newline & & & \ddots & \ddots & \ddots & \newline & & & & 1 & -2 & 1 \newline & & & & & 1 & -2 \newline \end{array} \right ) \end{equation}\begin{equation} x = \left ( \begin{array}{c} f_\mathrm{1} \\\ f_\mathrm{2} \\\ f_\mathrm{3} \\\ \vdots \\\ \vdots \\\ f_\mathrm{N-2} \\\ f_\mathrm{N-1} \\\ \end{array} \right ) \end{equation}\begin{equation} b = \Delta x^2 \left ( \begin{array}{c} g_\mathrm{1} \\\ g_\mathrm{2} \\\ g_\mathrm{3} \\\ \vdots \\\ \vdots \\\ g_\mathrm{N-2} \\\ g_\mathrm{N-1}\\\ \end{array} \right ) \end{equation}

Then we just solve $A x = b$

In [13]:

import scipy.linalg as linalg

# our grid -- including endpoints
N = 100
x = np.linspace(0.0, 2.0*np.pi, N, endpoint=True)
dx = x[1]-x[0]

# our source
g = np.sin(x)

# our matrix will be tridiagonal, with [1, -2, 1] on the diagonals
# we only solve for the N-2 interior points

# diagonal
d = -2*np.ones(N-2)

# upper -- note that the upper diagonal has 1 less element than the
# main diagonal.  The SciPy banded solver wants the matrix in the 
# form:
#
# *    a01  a12  a23  a34  a45    <- upper diagonal
# a00  a11  a22  a33  a44  a55    <- diagonal
# a10  a21  a32  a43  a54   *     <- lower diagonal
#

u = np.ones(N-2)
u[0] = 0.0

# lower
l = np.ones(N-2)
l[N-3] = 0.0

# put the upper, diagonal, and lower parts together as a banded matrix
A = np.matrix([u,d,l])

# solve A sol = dx**2 g for the inner N-2 points
sol = linalg.solve_banded((1,1), A, dx**2*g[1:N-1])

pylab.plot(x[1:N-1], sol)

Out[13]:

[<matplotlib.lines.Line2D at 0x1f21a2cf4f0>]

Eigenvalues¶

We can use np.linalg.eig to compute eigenvalues and eigenvectors of a matrix. It returns the eigenvalues and eigenvectors.

In [14]:

A = np.array([[1, -1], [1, 1]])
w, v = np.linalg.eig(A)
print(w)
print(v)

[1.+1.j 1.-1.j]
[[0.70710678+0.j         0.70710678-0.j        ]
 [0.        -0.70710678j 0.        +0.70710678j]]

Here, each entry in w is an eigenvalue of the matrix A, and each column vector in v is its corresponding eigenvector (if exists), such that a @ v[:,i] = w[i] * v[:,i].

In [15]:

A @ v[:,1] - w[1] * v[:,1]

Out[15]:

array([0.+0.j, 0.+0.j])

In [16]:

np.linalg.eig?

Signature: np.linalg.eig(a)
Docstring:
Compute the eigenvalues and right eigenvectors of a square array.

Parameters
----------
a : (..., M, M) array
    Matrices for which the eigenvalues and right eigenvectors will
    be computed

Returns
-------
w : (..., M) array
    The eigenvalues, each repeated according to its multiplicity.
    The eigenvalues are not necessarily ordered. The resulting
    array will be of complex type, unless the imaginary part is
    zero in which case it will be cast to a real type. When `a`
    is real the resulting eigenvalues will be real (0 imaginary
    part) or occur in conjugate pairs

v : (..., M, M) array
    The normalized (unit "length") eigenvectors, such that the
    column ``v[:,i]`` is the eigenvector corresponding to the
    eigenvalue ``w[i]``.

Raises
------
LinAlgError
    If the eigenvalue computation does not converge.

See Also
--------
eigvals : eigenvalues of a non-symmetric array.
eigh : eigenvalues and eigenvectors of a real symmetric or complex
       Hermitian (conjugate symmetric) array.
eigvalsh : eigenvalues of a real symmetric or complex Hermitian
           (conjugate symmetric) array.
scipy.linalg.eig : Similar function in SciPy that also solves the
                   generalized eigenvalue problem.
scipy.linalg.schur : Best choice for unitary and other non-Hermitian
                     normal matrices.

Notes
-----

.. versionadded:: 1.8.0

Broadcasting rules apply, see the `numpy.linalg` documentation for
details.

This is implemented using the ``_geev`` LAPACK routines which compute
the eigenvalues and eigenvectors of general square arrays.

The number `w` is an eigenvalue of `a` if there exists a vector
`v` such that ``a @ v = w * v``. Thus, the arrays `a`, `w`, and
`v` satisfy the equations ``a @ v[:,i] = w[i] * v[:,i]``
for :math:`i \in \{0,...,M-1\}`.

The array `v` of eigenvectors may not be of maximum rank, that is, some
of the columns may be linearly dependent, although round-off error may
obscure that fact. If the eigenvalues are all different, then theoretically
the eigenvectors are linearly independent and `a` can be diagonalized by
a similarity transformation using `v`, i.e, ``inv(v) @ a @ v`` is diagonal.

For non-Hermitian normal matrices the SciPy function `scipy.linalg.schur`
is preferred because the matrix `v` is guaranteed to be unitary, which is
not the case when using `eig`. The Schur factorization produces an
upper triangular matrix rather than a diagonal matrix, but for normal
matrices only the diagonal of the upper triangular matrix is needed, the
rest is roundoff error.

Finally, it is emphasized that `v` consists of the *right* (as in
right-hand side) eigenvectors of `a`.  A vector `y` satisfying
``y.T @ a = z * y.T`` for some number `z` is called a *left*
eigenvector of `a`, and, in general, the left and right eigenvectors
of a matrix are not necessarily the (perhaps conjugate) transposes
of each other.

References
----------
G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando, FL,
Academic Press, Inc., 1980, Various pp.

Examples
--------
>>> from numpy import linalg as LA

(Almost) trivial example with real e-values and e-vectors.

>>> w, v = LA.eig(np.diag((1, 2, 3)))
>>> w; v
array([1., 2., 3.])
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

Real matrix possessing complex e-values and e-vectors; note that the
e-values are complex conjugates of each other.

>>> w, v = LA.eig(np.array([[1, -1], [1, 1]]))
>>> w; v
array([1.+1.j, 1.-1.j])
array([[0.70710678+0.j        , 0.70710678-0.j        ],
       [0.        -0.70710678j, 0.        +0.70710678j]])

Complex-valued matrix with real e-values (but complex-valued e-vectors);
note that ``a.conj().T == a``, i.e., `a` is Hermitian.

>>> a = np.array([[1, 1j], [-1j, 1]])
>>> w, v = LA.eig(a)
>>> w; v
array([2.+0.j, 0.+0.j])
array([[ 0.        +0.70710678j,  0.70710678+0.j        ], # may vary
       [ 0.70710678+0.j        , -0.        +0.70710678j]])

Be careful about round-off error!

>>> a = np.array([[1 + 1e-9, 0], [0, 1 - 1e-9]])
>>> # Theor. e-values are 1 +/- 1e-9
>>> w, v = LA.eig(a)
>>> w; v
array([1., 1.])
array([[1., 0.],
       [0., 1.]])
File:      c:\users\win10hx64\appdata\local\programs\python\python310\lib\site-packages\numpy\linalg\linalg.py
Type:      function

In [ ]: