Linear algebra in Sympy/Python

Posted by ECON爱好者 on January 23, 2017   linear algbra   sympy
linearalgebra

The most important intuitioin of Linear algebra

https://youtu.be/kYB8IZa5AuE?t=74

In [ ]:
In [28]:
import matplotlib.pyplot as plt
%matplotlib inline
In [13]:
from sympy import *
In [3]:
from sympy import init_printing
init_printing()
In [2]:
from sympy import Matrix
In [ ]:
In [ ]:
In [48]:
x =  Matrix([1,3])
x
Out[48]:
$$\left[\begin{matrix}1\\3\end{matrix}\right]$$
In [49]:
2*x
Out[49]:
$$\left[\begin{matrix}2\\6\end{matrix}\right]$$
In [46]:
#fig, ax = plt.subplots()
fig = plt.figure()
ax = fig.add_subplot(111)
# Set the axes through the origin
for spine in ['left', 'bottom']:
    ax.spines[spine].set_position('zero')
for spine in ['right', 'top']:
    ax.spines[spine].set_color('none')


ax.set_xlim(-5, 5)
ax.set_ylim(-5, 7)
ax.grid()
vecs = ( (2, 6) , (1, 3))
for v in vecs:
    ax.annotate('', xy=v, xytext=(0, 0),
                arrowprops=dict(facecolor='red',
                shrink=0,
                alpha=0.7,
                width=0.5))
    ax.text(1.1 * v[0], 1.1 * v[1], str(v))
plt.show()
In [36]:
A = Matrix( [[ 2,1],
[-1, 1]] )

A
Out[36]:
$$\left[\begin{matrix}2 & 1\\-1 & 1\end{matrix}\right]$$
In [39]:
y = Matrix([5,2])
y
Out[39]:
$$\left[\begin{matrix}5\\2\end{matrix}\right]$$
In [40]:
A*x
Out[40]:
$$\left[\begin{matrix}5\\2\end{matrix}\right]$$
In [34]:
#fig, ax = plt.subplots()
fig = plt.figure()
ax = fig.add_subplot(111)
# Set the axes through the origin
for spine in ['left', 'bottom']:
    ax.spines[spine].set_position('zero')
for spine in ['right', 'top']:
    ax.spines[spine].set_color('none')


ax.set_xlim(-5, 5)
ax.set_ylim(-5, 5)
ax.grid()
vecs = ((1, 3), (5, 2))
for v in vecs:
    ax.annotate('', xy=v, xytext=(0, 0),
                arrowprops=dict(facecolor='blue',
                shrink=0,
                alpha=0.7,
                width=0.5))
    ax.text(1.1 * v[0], 1.1 * v[1], str(v))
plt.show()
In [51]:
phi=symbols('phi')
rotation=Matrix([[cos(phi), -sin(phi)],
                 [sin(phi), cos(phi)]])

rotation
Out[51]:
$$\left[\begin{matrix}\cos{\left (\phi \right )} & - \sin{\left (\phi \right )}\\\sin{\left (\phi \right )} & \cos{\left (\phi \right )}\end{matrix}\right]$$
Rotate
$$\begin{bmatrix} cos(\theta) & -sin(\theta)\\ sin(\theta) & cos(\theta) \\ \end{bmatrix}$$
stretching
$$\begin{bmatrix} \alpha & 0 \\ 0 & \alpha \\ \end{bmatrix}$$
In [ ]:
In [ ]:
In [41]:
#fig, ax = plt.subplots()
fig = plt.figure()
ax = fig.add_subplot(111)
# Set the axes through the origin
for spine in ['left', 'bottom']:
    ax.spines[spine].set_position('zero')
for spine in ['right', 'top']:
    ax.spines[spine].set_color('none')


ax.set_xlim(-5, 5)
ax.set_ylim(-5, 5)
ax.grid()
vecs = ((1, 0), (0, 1), (2,-1 ),(1,1))
for v in vecs:
    ax.annotate('', xy=v, xytext=(0, 0),
                arrowprops=dict(facecolor='green',
                shrink=0,
                alpha=0.7,
                width=0.5))
    ax.text(1.1 * v[0], 1.1 * v[1], str(v))
plt.show()
In [ ]:
In [62]:
M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
v = [1, 1, 1]
In [63]:
M
Out[63]:
$$\left[\begin{matrix}1 & 2 & 3\\4 & 5 & 6\\7 & 8 & 9\end{matrix}\right]$$
In [64]:
v
Out[64]:
$$\left [ 1, \quad 1, \quad 1\right ]$$
In [65]:
M.row(0).dot(v)
Out[65]:
$$6$$
In [66]:
M.col(0).dot(v)
Out[66]:
$$12$$
In [67]:
M.dot(v)
Out[67]:
$$\left [ 6, \quad 15, \quad 24\right ]$$
In [ ]:
In [2]:
A = Matrix( [[ 2,-3,-8, 7],
[-2,-1, 2,-7],
[ 1, 0,-3, 6]] )
In [3]:
A[0,1] # row 0, col 1of A
Out[3]:
-3
In [7]:
A[0:2,0:3] # top-left 2x3 submatrix of A
Out[7]:
$$\left[\begin{matrix}2 & -3 & -8\\-2 & -1 & 2\end{matrix}\right]$$
In [11]:
eye(2) # 2x2 identity matrix
Out[11]:
$$\left[\begin{matrix}1 & 0\\0 & 1\end{matrix}\right]$$
In [16]:
zeros(3,2)
Out[16]:
$$\left[\begin{matrix}0 & 0\\0 & 0\\0 & 0\end{matrix}\right]$$
In [17]:
A.transpose() # the same as A.T
Out[17]:
$$\left[\begin{matrix}2 & -2 & 1\\-3 & -1 & 0\\-8 & 2 & -3\\7 & -7 & 6\end{matrix}\right]$$

row operation

In [18]:
M = eye(3)
M[1,:] = M[1,:] + 3*M[0,:]
M
Out[18]:
$$\left[\begin{matrix}1 & 0 & 0\\3 & 1 & 0\\0 & 0 & 1\end{matrix}\right]$$
In [ ]:
In [26]:
import numpy as np
In [55]:
a=np.array([[1,2],[3,4]])
b=np.array([[11,12],[13,14]])
np.dot(a,b)
Out[55]:
array([[37, 40],
       [85, 92]])
In [56]:
np.inner(a,b)
Out[56]:
array([[35, 41],
       [81, 95]])
In [57]:
np.outer(a,b)
Out[57]:
array([[11, 12, 13, 14],
       [22, 24, 26, 28],
       [33, 36, 39, 42],
       [44, 48, 52, 56]])
In [ ]:
In [ ]:
In [29]:
A = Matrix(a)
A
Out[29]:
$$\left[\begin{matrix}1 & 2\\3 & 4\end{matrix}\right]$$
In [31]:
B = Matrix(b)
B
Out[31]:
$$\left[\begin{matrix}11 & 12\\13 & 14\end{matrix}\right]$$
In [32]:
A*B # dot
Out[32]:
$$\left[\begin{matrix}37 & 40\\85 & 92\end{matrix}\right]$$
In [ ]:
In [ ]:
In [61]:
A = Matrix([ [1, 3, 4],
[2, 9, 14],
[5, 12, 18]])
A
Out[61]:
$$\left[\begin{matrix}1 & 3 & 4\\2 & 9 & 14\\5 & 12 & 18\end{matrix}\right]$$
In [62]:
b = Matrix([8, 25,39])
b.T
Out[62]:
$$\left[\begin{matrix}8 & 25 & 39\end{matrix}\right]$$
In [63]:
b
Out[63]:
$$\left[\begin{matrix}8\\25\\39\end{matrix}\right]$$
In [65]:
A.row_join(b)
Out[65]:
$$\left[\begin{matrix}1 & 3 & 4 & 8\\2 & 9 & 14 & 25\\5 & 12 & 18 & 39\end{matrix}\right]$$
In [66]:
A.row_join(b).rref()
Out[66]:
$$\left ( \left[\begin{matrix}1 & 0 & 0 & 3\\0 & 1 & 0 & -1\\0 & 0 & 1 & 2\end{matrix}\right], \quad \left [ 0, \quad 1, \quad 2\right ]\right )$$
In [ ]:
In [71]:
# Gaussian_elimination
A = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 10]])
b = Matrix([3, 6, 9])
sol, params = A.gauss_jordan_solve(b)
sol
Out[71]:
$$\left[\begin{matrix}-1\\2\\0\end{matrix}\right]$$

Determinant

http://www.mathsisfun.com/algebra/matrix-determinant.html

In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [59]:
A = Matrix([ [1, 1, -1],
[2, 1, 1],
[2, 3, -5]])
A
Out[59]:
$$\left[\begin{matrix}1 & 1 & -1\\2 & 1 & 1\\2 & 3 & -5\end{matrix}\right]$$
In [61]:
A.rref()
Out[61]:
$$\left ( \left[\begin{matrix}1 & 0 & 2\\0 & 1 & -3\\0 & 0 & 0\end{matrix}\right], \quad \left [ 0, \quad 1\right ]\right )$$
In [60]:
A.rank()
Out[60]:
$$2$$
In [ ]:
In [68]:
b = Matrix([-3, -1,-10])
b
Out[68]:
$$\left[\begin{matrix}-3\\-1\\-10\end{matrix}\right]$$
In [69]:
A.row_join(b)
Out[69]:
$$\left[\begin{matrix}1 & 1 & -1 & -3\\2 & 1 & 1 & -1\\2 & 3 & -5 & -10\end{matrix}\right]$$
In [70]:
A.row_join(b).rref()
Out[70]:
$$\left ( \left[\begin{matrix}1 & 0 & 2 & 0\\0 & 1 & -3 & 0\\0 & 0 & 0 & 1\end{matrix}\right], \quad \left [ 0, \quad 1, \quad 3\right ]\right )$$
In [ ]:
In [ ]:
In [75]:
M = Matrix( [[1, 2, 3],
[2,-2, 4],
[2, 2, 5]] )
M
Out[75]:
$$\left[\begin{matrix}1 & 2 & 3\\2 & -2 & 4\\2 & 2 & 5\end{matrix}\right]$$
In [76]:
M.det()
Out[76]:
$$2$$
In [72]:
A = Matrix( [[1,2],
[3,9]] )
A
Out[72]:
$$\left[\begin{matrix}1 & 2\\3 & 9\end{matrix}\right]$$
In [73]:
A.inv()
Out[73]:
$$\left[\begin{matrix}3 & - \frac{2}{3}\\-1 & \frac{1}{3}\end{matrix}\right]$$
In [74]:
A.inv()*A
Out[74]:
$$\left[\begin{matrix}1 & 0\\0 & 1\end{matrix}\right]$$
In [9]:
A*A.inv()
Out[9]:
$$\left[\begin{matrix}1 & 0\\0 & 1\end{matrix}\right]$$
In [ ]:
In [57]:
Out[57]:
array([[ 1.,  1.,  1.],
       [ 1.,  1.,  1.],
       [ 1.,  1.,  1.],
       [ 1.,  1.,  1.]])
In [58]:
Out[58]:
$$2$$

What is an Eigenvector? Eigenvalue?

https://www.youtube.com/watch?v=wXCRcnbCsJA

Eigenvectors and eigenvalues | Essence of linear algebra, chapter 10

https://www.youtube.com/watch?v=PFDu9oVAE-g&t=683s

In [17]:
A = Matrix( [[ 9, -2],
[-2, 6]] )
A.eigenvals() # same as solve( det(A-eye(2)*x), x)
Out[17]:
$$\left \{ 5 : 1, \quad 10 : 1\right \}$$
In [ ]:
A.eigenvects()
In [ ]:
In [15]:
x = symbols('x')
In [ ]:
In [16]:
solve( det(A-eye(2)*x), x)
Out[16]:
$$\left [ 5, \quad 10\right ]$$
In [ ]:
In [79]:
import matplotlib.pyplot as plt
import numpy as np
from scipy.linalg import eig

A = ((1, 2),
     (2, 1))
A = np.array(A)
evals, evecs = eig(A)
evecs = evecs[:, 0], evecs[:, 1]
In [80]:
A
Out[80]:
array([[1, 2],
       [2, 1]])
In [81]:
evecs
Out[81]:
(array([ 0.70710678,  0.70710678]), array([-0.70710678,  0.70710678]))
In [82]:
evals
Out[82]:
array([ 3.+0.j, -1.+0.j])
In [ ]:
In [77]:
fig, ax = plt.subplots()
# Set the axes through the origin
for spine in ['left', 'bottom']:
    ax.spines[spine].set_position('zero')
for spine in ['right', 'top']:
    ax.spines[spine].set_color('none')
ax.grid(alpha=0.4)

xmin, xmax = -3, 3
ymin, ymax = -3, 3
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
# ax.set_xticks(())
# ax.set_yticks(())

# Plot each eigenvector
for v in evecs:
    ax.annotate('', xy=v, xytext=(0, 0),
                arrowprops=dict(facecolor='blue',
                shrink=0,
                alpha=0.6,
                width=0.5))

# Plot the image of each eigenvector
for v in evecs:
    v = np.dot(A, v)
    ax.annotate('', xy=v, xytext=(0, 0),
                arrowprops=dict(facecolor='red',
                shrink=0,
                alpha=0.6,
                width=0.5))

# Plot the lines they run through
x = np.linspace(xmin, xmax, 3)
for v in evecs:
    a = v[1] / v[0]
    ax.plot(x, a * x, 'b-', lw=0.4)


plt.show()
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [18]:
Q, L = A.diagonalize()
Q # the matrix of eigenvectors
Out[18]:
$$\left[\begin{matrix}1 & -2\\2 & 1\end{matrix}\right]$$
In [19]:
Q.inv()
Out[19]:
$$\left[\begin{matrix}\frac{1}{5} & \frac{2}{5}\\- \frac{2}{5} & \frac{1}{5}\end{matrix}\right]$$
In [20]:
Q*L*Q.inv() # eigendecomposition of A
Out[20]:
$$\left[\begin{matrix}9 & -2\\-2 & 6\end{matrix}\right]$$
In [ ]:
In [10]:
# http://emweb.unl.edu/Math/mathweb/vectors/vectors.html
v1 = Matrix([1,2,3])
v2 = Matrix([4,5,6])
v3 = v1.cross(v2)
v3
Out[10]:
$$\left[\begin{matrix}-3\\6\\-3\end{matrix}\right]$$
In [40]:
#A.dot(B)
Out[40]:
$$\left [ 47, \quad 55, \quad 70, \quad 82\right ]$$
In [ ]:
Reduced row echelon form

The Gauss–Jordan elimination procedure is a sequence of row operations you can perform on any matrix to bring it to its reduced row echelon form (RREF). In SymPy, matrices have a rref method that computes their RREF:

In [22]:
A = Matrix( [[2,-3,-8, 7],
[-2,-1,2,-7],
[1 ,0,-3, 6]])
A
Out[22]:
$$\left[\begin{matrix}2 & -3 & -8 & 7\\-2 & -1 & 2 & -7\\1 & 0 & -3 & 6\end{matrix}\right]$$
In [23]:
A.rref()
Out[23]:
$$\left ( \left[\begin{matrix}1 & 0 & 0 & 0\\0 & 1 & 0 & 3\\0 & 0 & 1 & -2\end{matrix}\right], \quad \left [ 0, \quad 1, \quad 2\right ]\right )$$

Note the rref method returns a tuple of values: the first value is the RREF of A, while the second tells you the indices of the leading ones (also known as pivots) in the RREF of A. To get just the RREF of A, select the 0th entry form the tuple: A.rref()[0].

In [ ]:

sympy: Symbolic Mathematics in Python

http://byumcl.bitbucket.org/bootcamp2013/labs/sympy.html

Scientific Computing with Python 3

http://www.bogotobogo.com/python/python_numpy_matrix_tutorial.php

In [42]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
from scipy.interpolate import interp2d

fig = plt.figure()
ax = fig.gca(projection='3d')

x_min, x_max = -5, 5
y_min, y_max = -5, 5

alpha, beta = 0.2, 0.1

ax.set_xlim((x_min, x_max))
ax.set_ylim((x_min, x_max))
ax.set_zlim((x_min, x_max))

# Axes
ax.set_xticks((0,))
ax.set_yticks((0,))
ax.set_zticks((0,))
gs = 3
z = np.linspace(x_min, x_max, gs)
x = np.zeros(gs)
y = np.zeros(gs)
ax.plot(x, y, z, 'k-', lw=2, alpha=0.5)
ax.plot(z, x, y, 'k-', lw=2, alpha=0.5)
ax.plot(y, z, x, 'k-', lw=2, alpha=0.5)


# Fixed linear function, to generate a plane
def f(x, y):
    return alpha * x + beta * y

# Vector locations, by coordinate
x_coords = np.array((3, 3))
y_coords = np.array((4, -4))
z = f(x_coords, y_coords)
for i in (0, 1):
    ax.text(x_coords[i], y_coords[i], z[i], r'$a_{}$'.format(i+1), fontsize=14)

# Lines to vectors
for i in (0, 1):
    x = (0, x_coords[i])
    y = (0, y_coords[i])
    z = (0, f(x_coords[i], y_coords[i]))
    ax.plot(x, y, z, 'b-', lw=1.5, alpha=0.6)


# Draw the plane
grid_size = 20
xr2 = np.linspace(x_min, x_max, grid_size)
yr2 = np.linspace(y_min, y_max, grid_size)
x2, y2 = np.meshgrid(xr2, yr2)
z2 = f(x2, y2)
ax.plot_surface(x2, y2, z2, rstride=1, cstride=1, cmap=cm.jet,
                linewidth=0, antialiased=True, alpha=0.2)
plt.show()
In [54]:
x = symbols('x')
plot((x**2, (x, -6, 6)), (x, (x, -5, 5)))
Out[54]:
<sympy.plotting.plot.Plot at 0x7f9c86e62cd0>
In [55]:
Out[55]:
<sympy.plotting.plot.Plot at 0x7f9c85e333d0>
In [56]:
In [ ]:

Related Posts