Python (SymPy) Programming for matrix calculations

Here I'd like to share how to deal with matrix calculation with Python (SymPy). For an introduction to how to use SymPy, see

pianofisica.hatenablog.com

Matri manipulation

Matri manipulation

Input matrices

To input matrices in Python (SymPy), put the following

import sympy
sympy.init_printing()
sympy.var('m11, m12, n11, n21')
sympy.var('a11, a12, a21, a22, b11, b12, b21, b22') 
M = sympy.Matrix([
[m11,m12]
])
N = sympy.Matrix([
[n11],
[n21]
])
A = sympy.Matrix([
[a11,a12],
[a21,a22]
])
B = sympy.Matrix([
[b11,b12],
[b21,b22]
])

Here $M$ is a matrix of type- $(1,2)$ , $N$ is of $(2,1)$ , $A$ and $B$ are of $(2,2)$ :

$\displaystyle{M=\begin{bmatrix} m_{11} & m_{12} \end{bmatrix} \qquad N=\begin{bmatrix} n_{11} \\ n_{21} \end{bmatrix} \qquad A=\begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix} \qquad B=\begin{bmatrix} b_{11} & b_{12} \\ b_{21} & b_{22} \end{bmatrix}}$

Refer matrix elements

A[0,1]   #  (k ,l) component of matrix A is given by  A[ k-1 , l-1 ]

B[1,1]

and so forth.

Operations of matrices (Product, Sum, Scalar multiplication, Power)

By definition of matrix product, for the case of our settings, matrix products
$AB$ , $MA$ , $BN$ , $MBN$ , and etcetera are allowed:

A*B

M*A

B*N

M*B*N

Power of matrices is given by

A**2

A linear combination of the matrices of the same type is obtained by

sympy.var('a, b')
a*A+b*B

Find inverse matrix

A.inv()

A**(-1)

Solve simultaneous linear equations

As an application of inverse matrix, let us consider simultaneous linear equations
$~ a_{11}x_1 + a_{12} x_2 + \dots + a_{1n} x_n=y_1 \\ ~ a_{21}x_1 + a_{22} x_2 + \dots + a_{2n} x_n=y_2 \\ ~~ \vdots \\ ~ a_{n1}x_1 + a_{n2} x_2 + \dots + a_{nn} x_n=y_n$

This is represented in terms of matrices as

$\displaystyle{~\begin{pmatrix} a_{11} & a_{12} & \dots & a_{1n} \\ a_{21} & a_{22} & \dots & a_{2n} \\ \vdots & \vdots & \dots & \vdots \\ a_{n1} & a_{n2} & \dots & a_{nn} \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \\ \dots \\ x_n \end{pmatrix}=\begin{pmatrix} y_1 \\ y_2 \\ \dots \\ y_n \end{pmatrix}}$

Its solution is given by an inverse matrix of

$\displaystyle{~A_n=\begin{pmatrix} a_{11} & a_{12} & \dots & a_{1n} \\ a_{21} & a_{22} & \dots & a_{2n} \\ \vdots & \vdots & \dots & \vdots \\ a_{n1} & a_{n2} & \dots & a_{nn} \end{pmatrix}}$

(here we assume the existence of $A_n^{-1}$ for simplicity) by applying which, one finds

$\displaystyle{~\begin{pmatrix} x_1 \\ x_2 \\ \dots \\ x_n \end{pmatrix}=A_n^{-1}\begin{pmatrix} y_1 \\ y_2 \\ \dots \\ y_n \end{pmatrix}}$

For example, simultaneous equations for three variables $x$ , $y$ and $z$

$~ ax+by+cz=j \\ ~ dx+ey+fz=k \\ ~ gx+hy+iz=l$

is solved by finding the inverse of the matrix

$\displaystyle{~A_3=\begin{pmatrix} a & b & c \\ d & e & f \\ g & h & i \end{pmatrix}}$

import sympy
sympy.init_printing()
sympy.var('a, b, c, d, e, f, g, h, i') 
A3 = sympy.Matrix([
[a, b, c],
[d, e, f],
[g, h, i]
])
sympy.cancel(A3**(-1))

then the inverse is found to be

$\displaystyle{~A_3^{-1}=\frac{1}{\det A_3}\begin{pmatrix} ei-fh & -(bi-ch) & bf-ce \\ -(𝑑𝑖−𝑓𝑔) & 𝑎𝑖−𝑐𝑔 & -(𝑎𝑓−𝑐𝑑) \\ 𝑑ℎ−𝑒𝑔 & -(𝑎ℎ−𝑏𝑔) & 𝑎𝑒−𝑏𝑑 \end{pmatrix}}$

where

$\displaystyle{~\det A_3=a(ei-fh)-b(di-fg)+c(dh-eg)}$

Finally, take its product with the vector

$\displaystyle{~Y=\begin{pmatrix} j \\ k \\ l \end{pmatrix}}$

one finds the solution

$\displaystyle{~\begin{pmatrix} x \\ y \\ z \end{pmatrix} = A_3^{-1}Y= \frac{1}{\det A_3}\begin{pmatrix} j(ei-fh) -k(bi-ch) +l(bf-ce) \\ -j(𝑑𝑖−𝑓𝑔) +k( 𝑎𝑖−𝑐𝑔) -l(𝑎𝑓−𝑐𝑑) \\ j(𝑑ℎ−𝑒𝑔) -k(𝑎ℎ−𝑏𝑔) +l( 𝑎𝑒−𝑏𝑑) \end{pmatrix}}$

To make a comparison with the result using "solve", put

eq1=sympy.Eq(a*x+b*y+c*z, j)
eq2=sympy.Eq(d*x+e*y+f*z, k)
eq3=sympy.Eq(g*x+h*y+i*z, l)
sympy.solve ([eq1, eq2, eq3], [x, y, z])

then

$\displaystyle{x={{j(ei-fh)-k(bi-ch)+l(bf-ce)}\over{a(ei-fh)-b(di-fg)+c(dh-eg)}} \\ y={{−𝑗(𝑑𝑖−𝑓𝑔)+𝑘(𝑎𝑖−𝑐𝑔)−𝑙(𝑎𝑓−𝑐𝑑)}\over{a(ei-fh)-b(di-fg)+c(dh-eg)}} \\ z={{𝑗(𝑑ℎ−𝑒𝑔)−𝑘(𝑎ℎ−𝑏𝑔)+𝑙(𝑎𝑒−𝑏𝑑)}\over{a(ei-fh)-b(di-fg)+c(dh-eg)} }}$

and thus they coincide.

Determinant

Determinant of matrix $A$ is calculated by

A.det()

$\widetilde{A}$ , the minor of $A$ is

A.adjugate()

indeed from

A.inv() == 1/(A.det())*A.adjugate()

thus

$\displaystyle{\widetilde{A}=\left(\det A\right) A^{-1}}$

Identity matrix

Identity matrix is

sympy.eye(2)

The argument specifies the matrix size.

For example,

sympy.simplify(A.inv()*A) == sympy.eye(2)

Trace

A.trace()

Transposed matrix

The transposed matrix of $A$ defined by

$\displaystyle{A^t=\begin{bmatrix} a_{11} & a_{21} \\ a_{12} & a_{22} \end{bmatrix}}$

is obtained by

A.transpose()

Adjoint matrix

The adjoint matrix $A$

$\displaystyle{A^\dagger=\begin{bmatrix} \overline{a_{11}} & \overline{a_{21}} \\ \overline{a_{12}} & \overline{a_{22}} \end{bmatrix}}$

is done by

A.adjoint()

Diagonalize matrices

Find eigenvalues

A.eigenvals()

Find eigenvectors

A.eigenvects()

Diagonalized matrix

A.diagonalize()

Keywords：Python, SymPy, Matrix