Linear algebra examples

The following illustrates how to use the numpy and scipy.linalg modules to work with matricies by solving problem A.8 from Introduction to Quantum Mechanics, 2nd Edition by David Griffiths and another example eigenvalue problem.

In [1]:
import numpy as np
import scipy.linalg as la

Problem A.8 (Griffiths)

Given the matricies $\mathbf{A}$ and $\mathbf{B}$ compute the quantities specified in parts (a) through (i) below.

Define matrix $\mathbf{A}$ and $\mathbf{B}$.

In [2]:
A = np.array([[-1,1,1j],
              [2,0,3],
              [2j,-2j,2]])
B = np.array([[2,0,-1j],
              [0,1,0],
              [1j,3,2]])
print('A = ',A)
print('B = ',B)

A =  [[-1.+0.j  1.+0.j  0.+1.j]
 [ 2.+0.j  0.+0.j  3.+0.j]
 [ 0.+2.j -0.-2.j  2.+0.j]]
B =  [[ 2.+0.j  0.+0.j -0.-1.j]
 [ 0.+0.j  1.+0.j  0.+0.j]
 [ 0.+1.j  3.+0.j  2.+0.j]]

(a) Matrix addition A+B

In [3]:
C = A+B
print(C)

[[1.+0.j 1.+0.j 0.+0.j]
 [2.+0.j 1.+0.j 3.+0.j]
 [0.+3.j 3.-2.j 4.+0.j]]

(b) Matrix multiplication AB

The NumPy uses the @ symbol for the matrix multiplication operator.

In [4]:
A@B
Out[4]:
array([[-3.+0.j,  1.+3.j,  0.+3.j],
       [ 4.+3.j,  9.+0.j,  6.-2.j],
       [ 0.+6.j,  6.-2.j,  6.+0.j]])

(c) Commutator [A,B]

In [5]:
A@B-B@A
Out[5]:
array([[-3.+0.j,  1.+3.j,  0.+3.j],
       [ 2.+3.j,  9.+0.j,  3.-2.j],
       [-6.+3.j,  6.+1.j, -6.+0.j]])

(d) Transpose of A ($\tilde{\mathbf{A}}$)

In [6]:
A.T
Out[6]:
array([[-1.+0.j,  2.+0.j,  0.+2.j],
       [ 1.+0.j,  0.+0.j, -0.-2.j],
       [ 0.+1.j,  3.+0.j,  2.+0.j]])

(e) Complex conjugate of A ($\mathbf{A}^\ast$)

In [7]:
A.conjugate()
Out[7]:
array([[-1.-0.j,  1.-0.j,  0.-1.j],
       [ 2.-0.j,  0.-0.j,  3.-0.j],
       [ 0.-2.j, -0.+2.j,  2.-0.j]])

(f) Hermitian conjugate or adjoint of A ($\mathbf{A}^\dagger$)

In [8]:
A.conjugate().T
Out[8]:
array([[-1.-0.j,  2.-0.j,  0.-2.j],
       [ 1.-0.j,  0.-0.j, -0.+2.j],
       [ 0.-1.j,  3.-0.j,  2.-0.j]])

(g) Trace ($\mathrm{Tr}(\mathbf{B})$)

In [9]:
np.trace(B)
Out[9]:
(5+0j)

(h) Determinant of B ($\mathrm{det}(\mathbf{B})$)

In [10]:
la.det(B)
Out[10]:
(3+0j)

(i) Inverse of B ($\mathbf{B}^{-1}$)

Check that $\mathbf{BB}^{-1} = \mathbf{I}$. Does $\mathbf{A}$ have an inverse?

In [11]:
Binv = la.inv(B)
In [12]:
print(Binv)

[[ 0.66666667+0.j          0.        -1.j          0.        +0.33333333j]
 [ 0.        +0.j          1.        +0.j          0.        +0.j        ]
 [ 0.        -0.33333333j -2.        -0.j          0.66666667+0.j        ]]
In [13]:
BBI = Binv@B
print(BBI)

[[ 1.00000000e+00+0.00000000e+00j  0.00000000e+00-5.55111512e-17j
   0.00000000e+00+0.00000000e+00j]
 [ 0.00000000e+00+0.00000000e+00j  1.00000000e+00+0.00000000e+00j
   0.00000000e+00+0.00000000e+00j]
 [ 0.00000000e+00+0.00000000e+00j -1.11022302e-16+0.00000000e+00j
   1.00000000e+00+0.00000000e+00j]]

This confirms that $\mathbf{BB}^{-1} = \mathbf{I}$. The very small numbers in the off-diagonal elements are due to rounding errors.

If $\mathrm{det}(\mathbf{A}) \neq 0$, then $\mathbf{A}$ has an inverse.

In [14]:
la.det(A)
Out[14]:
0j

$\mathrm{det}(\mathbf{A}) = 0$ so $\mathbf{A}$ does not have an inverse.

Example eigenvalue and eigenvector problem

Find the eigenvalues and eigenvectors of the matrix $\mathbf{M}$ below.

In [15]:
M = np.array([[1,0,1],
              [0,1,0],
              [1,0,1]])

The la.eig() function takes a square matrix as its argument and returns two arrays. One is an a one-dimensional array that contains all the eigenvalues. The second array is a two-dimensional array with each column of the array being one of the eigenvectors. The code below illustrates how to extract the eigenvalues and eigenvectors from the output of la.eig().

In [16]:
e_val,e_vec = la.eig(M)

In this case e_val is an array with three elements—one element for each eigenvalue.

In [17]:
print('The eigenvalues are:',e_val)

The eigenvalues are: [2.+0.j 0.+0.j 1.+0.j]

As expected for this matrix, all of the eigenvalues are real. The imaginary part of each eigenvalue is zero.

Now we can extract the eigenvector for each eigenvalue and print the results.

In [18]:
for idx in range(len(e_val)):
    value = e_val[idx]
    vector = e_vec[:,idx]
    print(f'For the eigenvalue {value}, the eigenvector is {vector}.')

For the eigenvalue (2+0j), the eigenvector is [0.70710678 0.         0.70710678].
For the eigenvalue 0j, the eigenvector is [-0.70710678  0.          0.70710678].
For the eigenvalue (1+0j), the eigenvector is [0. 1. 0.].