--- author: - "[Eugeniy E. Mikhailov](http://physics.wm.edu/~evmik/) and [Greg Bentsen](https://physics.wm.edu/~gbentsen/)" title: "Numerical Optimization VI: Optimization and linear algebra" --- ## Today Linear algebra formulation - solving system of equation Bringing back optimization - solving over constrained system of equation ## Mobile problem ::: columns ::: {.column width="70% "} Someone provided us with 6 weights and 3 rods. We need to calculate the positions of suspension points to have a balanced system. If the system is in equilibrium, torque must be zero at every pivot point $$\begin{aligned} w_1 x_1 - (w_2 + w_3 + w_4 + w_5 + w_6) x_2 &=& 0 \\ w_3 x_3 - (w_4 + w_5 + w_6) x_4 &=& 0 \\ w_5 x_5 - w_6 x_6 &=& 0 \end{aligned}$$ Let's define $w_{26}=w_2 + w_3 + w_4 + w_5+w_6$ and $w_{46}=w_4+w_5+w_6$ $$\begin{aligned} w_1 x_1 - w_{26} x_2 &=& 0 \\ w_3 x_3 - w_{46} x_4 &=& 0 \\ w_5 x_5 - w_6 x_6 &=& 0 \\ \end{aligned}$$ We need 3 more equations. Let's constrain the length of every rod $$\begin{aligned} x_1 + x_2 &=& L_{12} \\ x_3 + x_4 &=& L_{34} \\ x_5 + x_6 &=& L_{56} \end{aligned}$$ ::: ::: {.column width="30%"} ![Mobile scketch](./pics/mobile.svg){width="100%"} ::: ::: ## Mobile problem continued ::: columns ::: {.column width="30%"} ![Mobile scketch](./pics/mobile.svg){width="100%"} ::: ::: {.column width="70% "} ::: columns ::: {.column width="35%"} $$\begin{aligned} w_1 x_1 - w_{26} x_2 &=& 0 \\ w_3 x_3 - w_{46} x_4 &=& 0 \\ w_5 x_5 - w_6 x_6 &=& 0 \\ x_1 + x_2 &=& L_{12} \\ x_3 + x_4 &=& L_{34} \\ x_5 + x_6 &=& L_{56} \end{aligned}$$ ::: ::: {.column width="55%"} $$\begin{aligned} \begin{pmatrix} w_1 & -w_{26} & 0 & 0 & 0 & 0 \\ 0 & 0 & w_3 & -w_{46} & 0 & 0 \\ 0 & 0 & 0 & 0 & w_5 & -w_6 \\ 1 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 1 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \\ x_5 \\ x_6 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ L_{12} \\ L_{34} \\ L_{56} \end{pmatrix} \end{aligned}$$ ::: definition $$\begin{aligned} \sum_j A_{ij} x_j = B_i \to {\mathbf A} {\mathbf x} = {\mathbf B} \end{aligned}$$ **Bold** marks matrices. Vector such as ${\mathbf x}$ and ${\mathbf B}$ are matrices too (1D). ::: ::: ::: Python has a lot of built-in functions to solve problems in this form ::: ::: ## Linear Algebra inverse matrix method ::: definition $$\begin{aligned} {\mathbf A} {\mathbf x} &=& {\mathbf B} \\ {\mathbf A^{-1}} {\mathbf A} {\mathbf x} &=& {\mathbf A^{-1}}{\mathbf B} \\ {\mathbf x} &=& {\mathbf A^{-1}}{\mathbf B} \end{aligned}$$ Analytical solution ${\mathbf x} = {\mathbf A^{-1}}{\mathbf B}$ exists only if $\text{determinant}({\mathbf A}) \ne 0$ Remember that ${\mathbf A^{-1}} {\mathbf A} = {\mathbf A} {\mathbf A^{-1}} = {\mathbf I}$, where ${\mathbf I}$ is the unity matrix. ::: ::: columns ::: {.column width="55%"} Python's straight forward implementation (not the fastest) ``` .python import numpy as np # define A and B marices invA = np.linalg.inv(A) # Pay attention @ is not the same as * # @ also is matrix multiplication by Linear Algebra rules # * is the element wise multiplication x = invA @ B ``` ::: ::: {.column width="45%"} Inverting matrix is computationally expensive! Python's faster way by using equation solver directly (recommended) ``` .python import numpy as np # define A and B marices invA = inv(A) x = np.linalg.solve(A, B) ``` ::: ::: ## Examples ::: columns ::: {.column width="55%"} $$\begin{aligned} 2 x_1 + 3 x_2 &= 8 \\ 5 x_1 + 1 x_2 &= 7 \end{aligned}$$ ::: ::: {.column width="55%"} $${\mathbf A} {\mathbf x} = {\mathbf B}$$ $$\begin{aligned} \begin{pmatrix} 2 & 3\\ 5 & 1 \end{pmatrix} \cdot \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} 8 \\ 7 \end{pmatrix} \end{aligned}$$ ::: ::: ::: columns ::: {.column width="55%"} ``` .python import numpy as np A = np.array([[ 2, 3], [5, 1]]) B = np.array([[8], [7]]) invA = np.linalg.inv(A) x = invA @ B In [9]: invA Out[9]: array([[-0.07692308, 0.23076923], [ 0.38461538, -0.15384615]]) In [10]: x Out[10]: array([[1.], [2.]]) ``` ::: ::: {.column width="45%"} ``` .python import numpy as np A = np.array([[ 2, 3], [5, 1]]) B = np.array([[8], [7]]) x = np.linalg.solve(A, B) In [12]: x Out[12]: array([[1.], [2.]]) # check In [15]: A @ x Out[15]: array([[8.], [7.]]) ``` ::: ::: ## Recall the mobile problem If $w_1=20$, $w_2=5$, $w_3=3$, $w_4=7$, $w_5=2$, $w_6=3$, $L_{12}=2$, $L_{34}=1$, and $L_{56}=3$, then $w_{26}=20$ and $w_{46}=12$. ::: columns ::: {.column width="35%"} $$\begin{aligned} \begin{pmatrix} 20 & -20 & 0 & 0 & 0 & 0 \\ 0 & 0 & 3 & -12 & 0 & 0 \\ 0 & 0 & 0 & 0 & 2 & -3 \\ 1 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 1 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \\ x_5 \\ x_6 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 2 \\ 1 \\ 3 \end{pmatrix} \end{aligned}$$ ::: ::: {.column width="35%"} ``` .python import numpy as np A=np.array([ [20, -20, 0, 0, 0, 0], [0, 0, 3, -12, 0, 0], [0, 0, 0, 0, 2, -3], [1, 1, 0, 0, 0, 0], [0, 0, 1, 1, 0, 0], [0, 0, 0, 0, 1, 1], ]) B= np.array([[0], [0], [0], [2], [1], [3] ]) x=np.linalg.inv(A) @ B ``` ::: ::: {.column width="25%"} ``` .python In [33]: x Out[33]: array([[1. ], [1. ], [0.8], [0.2], [1.8], [1.2]]) ``` Check ``` .python In [32]: A @ x - B Out[32]: array([[0.0000000e+00], [0.0000000e+00], [0.0000000e+00], [0.0000000e+00], [0.0000000e+00], [4.4408921e-16]]) ``` ::: ::: ## To do or not to do the inverse matrix calculation `\alert{Solutions based on the inverse matrix calculation}`{=tex} involve extra steps (unnecessary for solution) and, thus, `\alert{are slower}`{=tex} However, it is handy to calculate the inverse matrix in advance if you solve ${\mathbf A} {\mathbf x} = {\mathbf B}$ for different ${\mathbf B}$ with the same ${\mathbf A}$. ## Reformulating via optimization problem We can say that $${\mathbf A} {\mathbf x} = {\mathbf B}$$ is equivalent to ::: definition find ${\mathbf x}$ which minimizes $|{\mathbf A} {\mathbf x} - {\mathbf B}|$ ::: In python it can be done via ``` .python import numpy as np A = np.array([[ 2, 3], [5, 1]]) B = np.array([[8], [7]]) sol = np.linalg.lstsq(A, B) x = sol[0] In [40]: x Out[40]: array([[1.], [2.]]) ``` ## Optimization solver benefits Sounds useless but it allows to solve over constrained problems like $$\begin{aligned} 4x = 3 \\ 5x = 2 \end{aligned}$$ . . . Very common in physics. Mentally replace $x \to R$ $$\begin{aligned} I_1 R = V_1 \\ I_2 R = V_2 \\ \cdots \\ I_n R = V_n \\ \end{aligned}$$ . . . ``` .python import numpy as np A=np.array([ [4], [5] ]) B=np.array([ [3], [2] ]) sol = np.linalg.lstsq(A, B) x = sol[0] In [44]: x Out[44]: array([[0.53658537]]) ```