# Занятие 4
# Прикладная алгебра и численные методы
## Сингулярное разложение (SVD), линейная регрессия
https://numpy.org/doc/stable/reference/generated/numpy.linalg.svd.html#numpy.linalg.svd


```python
import numpy as np
import scipy.linalg
import sympy
import matplotlib.pyplot as plt
from copy import deepcopy
%matplotlib inline
```

## Сингулярное разложение (SVD)
$$
A = Q\Sigma P^*, \quad A_{m\times n},\ Q_{m\times m}, \ \Sigma_{m\times n}, \ P_{n\times n},
$$
$Q$, $P$ - ортогональные матрицы, $\Sigma$ - диагональная, на диагонали сингулярные числа.


## Пример 1
Найти SVD
$$
\left(
\begin{matrix}
1 & 0 & 0 & 1\\
0 & 1 & 0 & 1\\
0 & 0 & 1 & 1
\end{matrix}
\right)
$$
Вначале вычислим $A^*A$:


```python
A = sympy.Matrix([[1, 0, 0, 1], [0, 1, 0, 1], [0, 0, 1, 1]])
A_star_A = A.T*A
```

Получим собственные числа и собственные векторы с помощью eigenvects(), нормализуем векторы (чтобы норма была равна единице) методом normalized()


```python
A_star_A_sympy_ev = sympy.Matrix(A_star_A).eigenvects()
A_star_A_sympy_ev = [(item[0], item[1], [elem.normalized() for elem in item[2]]) for item in A_star_A_sympy_ev]
display(*A_star_A_sympy_ev)
```


    (0, 1, [Matrix([
      [-1/2],
      [-1/2],
      [-1/2],
      [ 1/2]])])



    (1, 2, [Matrix([
      [-sqrt(2)/2],
      [ sqrt(2)/2],
      [         0],
      [         0]]), Matrix([
      [-sqrt(2)/2],
      [         0],
      [ sqrt(2)/2],
      [         0]])])



    (4, 1, [Matrix([
      [sqrt(3)/6],
      [sqrt(3)/6],
      [sqrt(3)/6],
      [sqrt(3)/2]])])


Выделим собственные векторы, обозначим их e0, e11, e12, e4, они соответствуют собственным значеним 0, 1, 1 и 4. К двум векторам, соответствующим собственному значению 1 применим процесс ортогонализации Грамма-Шмидта 
$$
\begin{matrix}
e_1^{new} = e_1\\
e_2^{new} = e_2 - \frac{(e_1, e_2)}{(e_1, e_1)}e_1
\end{matrix}
$$
Полученный ортогональный вектор нормализуем, проверим ортогональность с помощью скалярного произведения:


```python
e0, e11, e12, e4 = A_star_A_sympy_ev[0][2] + A_star_A_sympy_ev[1][2] + A_star_A_sympy_ev[2][2]
e12 = (e12 - e11.dot(e12)*e11).normalized()
display(e11.dot(e12))
P = e4.row_join(e11).row_join(e12).row_join(e0)
display(P)
```


$\displaystyle 0$



$\displaystyle \left[\begin{matrix}\frac{\sqrt{3}}{6} & - \frac{\sqrt{2}}{2} & - \frac{\sqrt{6}}{6} & - \frac{1}{2}\\\frac{\sqrt{3}}{6} & \frac{\sqrt{2}}{2} & - \frac{\sqrt{6}}{6} & - \frac{1}{2}\\\frac{\sqrt{3}}{6} & 0 & \frac{\sqrt{6}}{3} & - \frac{1}{2}\\\frac{\sqrt{3}}{2} & 0 & 0 & \frac{1}{2}\end{matrix}\right]$


Построим векторы-столбцы матрицы $Q$ и проверим, что найдено разложение SVD для исходной матрицы:


```python
sigma = (2, 1, 1)
f1, f2, f3 = [A*ei/sigma[i] for i, ei in enumerate((e4, e11, e12))]
Q = f1.row_join(f2).row_join(f3)
Sig = sympy.Matrix([[2, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]])
display(Q, Sig, P, Q*Sig*P.T)
```


$\displaystyle \left[\begin{matrix}\frac{\sqrt{3}}{3} & - \frac{\sqrt{2}}{2} & - \frac{\sqrt{6}}{6}\\\frac{\sqrt{3}}{3} & \frac{\sqrt{2}}{2} & - \frac{\sqrt{6}}{6}\\\frac{\sqrt{3}}{3} & 0 & \frac{\sqrt{6}}{3}\end{matrix}\right]$



$\displaystyle \left[\begin{matrix}2 & 0 & 0 & 0\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\end{matrix}\right]$



$\displaystyle \left[\begin{matrix}\frac{\sqrt{3}}{6} & - \frac{\sqrt{2}}{2} & - \frac{\sqrt{6}}{6} & - \frac{1}{2}\\\frac{\sqrt{3}}{6} & \frac{\sqrt{2}}{2} & - \frac{\sqrt{6}}{6} & - \frac{1}{2}\\\frac{\sqrt{3}}{6} & 0 & \frac{\sqrt{6}}{3} & - \frac{1}{2}\\\frac{\sqrt{3}}{2} & 0 & 0 & \frac{1}{2}\end{matrix}\right]$



$\displaystyle \left[\begin{matrix}1 & 0 & 0 & 1\\0 & 1 & 0 & 1\\0 & 0 & 1 & 1\end{matrix}\right]$


Теперь то же самое, но с numpy, вычислим $A^*A$:


```python
A = np.array([[1, 0, 0, 1],
              [0, 1, 0, 1],
              [0, 0, 1, 1]])
A_star_A = np.matmul(A.T, A)
display(A_star_A)
```


    array([[1, 0, 0, 1],
           [0, 1, 0, 1],
           [0, 0, 1, 1],
           [1, 1, 1, 3]])


Найдем собственные числа и собственные векторы полученной матрицы:


```python
A_star_A_eigen_vals, A_star_A_eigen_vects = np.linalg.eig(A_star_A)
display('CЧ', A_star_A_eigen_vals, 'СВ', A_star_A_eigen_vects)
```


    'CЧ'



    array([ 4.000000e+00, -2.493665e-16,  1.000000e+00,  1.000000e+00])



    'СВ'



    array([[ 2.88675135e-01,  5.00000000e-01, -8.16496581e-01,
            -4.08248290e-01],
           [ 2.88675135e-01,  5.00000000e-01,  4.08248290e-01,
            -4.08248290e-01],
           [ 2.88675135e-01,  5.00000000e-01,  4.08248290e-01,
             8.16496581e-01],
           [ 8.66025404e-01, -5.00000000e-01,  3.34893759e-16,
             2.58096210e-16]])


Расположим сингуляные числа (квадратные корни из полученных собственных чисел) по убыванию, для этого сначала отсортируем их с помощью sort() по возрастанию, а затем запишем array в обратном порядке с помощью flip():


```python
A_star_A_eigen_vals.sort()
A_star_A_eigen_vals_reversed = np.flip(A_star_A_eigen_vals)
display(A_star_A_eigen_vals, A_star_A_eigen_vals_reversed)
```


    array([-2.493665e-16,  1.000000e+00,  1.000000e+00,  4.000000e+00])



    array([ 4.000000e+00,  1.000000e+00,  1.000000e+00, -2.493665e-16])


Обратите внимание, что .sort() изменяет array на месте, а flip() возвращает view записанного в обратном порядке array, не изменяя его.
## !!! 
По сути, мы получаем указатель на конец нашего array, так что все действия, которые мы проделаем с элементами A_star_A_eigen_vals_reversed автоматически распространятся на A_star_A_eigen_vals, поскольку это не два разных array, а один, только номера элементов считаются по-разному:


```python
arr1 = np.array([1, 2, 3, 4])
arr1_reversed = np.flip(arr1)
arr1[0] = 9
display(arr1, arr1_reversed)
arr1_reversed[-2] = 8
display(arr1, arr1_reversed)
```


    array([9, 2, 3, 4])



    array([4, 3, 2, 9])



    array([9, 8, 3, 4])



    array([4, 3, 8, 9])


Поскольку нам достаточно работать с A_star_A_eigen_vals_reversed, не будем делать deepcopy(), чтобы сохратить в неприкосновенности A_star_A_eigen_vals.

Осталось извлечь квадратные корни из положительных элементов A_star_A_eigen_vals_reversed, и получим невозрастающую последовательность сингулярных значений.


```python
sigmas = [round(np.sqrt(item), 1) for item in A_star_A_eigen_vals_reversed if item > 0]
sigmas
```




    [2.0, 1.0, 1.0]



Составим матрицу $\Sigma$:


```python
Sigma = np.hstack((np.diag(sigmas), np.zeros((3, 1))))
Sigma
```




    array([[2., 0., 0., 0.],
           [0., 1., 0., 0.],
           [0., 0., 1., 0.]])



Обратимся к полученным вместе с собственными числами собственным векторам:


```python
e4, e0, e11, e12 = [item.reshape((4, 1)) for item in  A_star_A_eigen_vects.T]
A_star_A_eigen_vects = [e4, e0, e11, e12]
display(A_star_A_eigen_vects)
```


    [array([[0.28867513],
            [0.28867513],
            [0.28867513],
            [0.8660254 ]]), array([[ 0.5],
            [ 0.5],
            [ 0.5],
            [-0.5]]), array([[-8.16496581e-01],
            [ 4.08248290e-01],
            [ 4.08248290e-01],
            [ 3.34893759e-16]]), array([[-4.08248290e-01],
            [-4.08248290e-01],
            [ 8.16496581e-01],
            [ 2.58096210e-16]])]


Вычислим нормы полученных векторов и скалярное произведение 


```python
display(*[np.linalg.norm(item) for item in A_star_A_eigen_vects], 'скалярное произведение', e11[0].dot(e12[0]))
```


    0.9999999999999999



    1.0



    1.0



    1.0



    'скалярное произведение'



    0.3333333333333335


Сначала заменим e12 на вектор, ортогональный e11, затем нормализуем векторы и составим из них матрицу $P$:


```python
A_star_A_eigen_vects = [item/np.linalg.norm(item) for item in A_star_A_eigen_vects]
e4, e0, e11, e12 = A_star_A_eigen_vects
e12 = e12 - (e11[0].dot(e12[0])/(e11[0].dot(e11[0])))*e11
print('(e11, e12new) =', e11[0].dot(e12[0]))
A_star_A_eigen_vects[-1] = e12
e4, e0, e11, e12 = [item/np.linalg.norm(item) for item in A_star_A_eigen_vects]
P = np.hstack((e4, e11, e12, e0))
P
```

    (e11, e12new) = -0.0





    array([[ 2.88675135e-01, -8.16496581e-01,  0.00000000e+00,
             5.00000000e-01],
           [ 2.88675135e-01,  4.08248290e-01, -7.07106781e-01,
             5.00000000e-01],
           [ 2.88675135e-01,  4.08248290e-01,  7.07106781e-01,
             5.00000000e-01],
           [ 8.66025404e-01,  3.34893759e-16,  1.04672831e-16,
            -5.00000000e-01]])



Составим матрицу $Q$:


```python
sigma = (2, 1, 1)
f1, f2, f3 = [np.matmul(A, ei)/sigma[i] for i, ei in enumerate((e4, e11, e12))]
Q = np.hstack((f1, f2, f3))
Sig = np.hstack((np.diag(sigma), np.zeros((3, 1))))
display(Q, Sig, P, np.matmul(np.matmul(Q, Sig), P.T))
```


    array([[ 5.77350269e-01, -8.16496581e-01,  1.04672831e-16],
           [ 5.77350269e-01,  4.08248290e-01, -7.07106781e-01],
           [ 5.77350269e-01,  4.08248290e-01,  7.07106781e-01]])



    array([[2., 0., 0., 0.],
           [0., 1., 0., 0.],
           [0., 0., 1., 0.]])



    array([[ 2.88675135e-01, -8.16496581e-01,  0.00000000e+00,
             5.00000000e-01],
           [ 2.88675135e-01,  4.08248290e-01, -7.07106781e-01,
             5.00000000e-01],
           [ 2.88675135e-01,  4.08248290e-01,  7.07106781e-01,
             5.00000000e-01],
           [ 8.66025404e-01,  3.34893759e-16,  1.04672831e-16,
            -5.00000000e-01]])



    array([[ 1.00000000e+00,  2.64089394e-16,  2.60350884e-16,
             1.00000000e+00],
           [ 9.67157293e-17,  1.00000000e+00, -1.88904597e-16,
             1.00000000e+00],
           [ 1.66992088e-16, -6.83580866e-17,  1.00000000e+00,
             1.00000000e+00]])


## Построение псевдообратной матрицы при помощи SVD
$$
A^+ = P\Sigma^+Q^*,\quad 
\Sigma^+ =
\left(
\begin{matrix}
\sigma_1^{-1} & ... & ... & ... & ... & 0\\
0 & \sigma_1^{-1} & ... & ... & ... & 0\\
0 & ... & ... & ... & ... & 0\\
0 & ... & ... & \sigma_r^{-1}  & ... & 0\\
0 & ... & ... & ... & ... & 0\\
\end{matrix}
\right)
$$


```python
Sigma_plus = np.vstack((np.diag([1/item for item in sigma]), np.zeros((1, 3))))
A_pinv_my = np.matmul(np.matmul(P, Sigma_plus), Q.transpose())
display(A_pinv_my, np.linalg.pinv(A))
```


    array([[ 0.75, -0.25, -0.25],
           [-0.25,  0.75, -0.25],
           [-0.25, -0.25,  0.75],
           [ 0.25,  0.25,  0.25]])



    array([[ 0.75, -0.25, -0.25],
           [-0.25,  0.75, -0.25],
           [-0.25, -0.25,  0.75],
           [ 0.25,  0.25,  0.25]])


## И наконец SVD от numpy:


```python
Q, sigma, P = np.linalg.svd(A, full_matrices=True)
Sig = np.hstack((np.diag(sigma), np.zeros((3, 1))))
display('P.T', P, 'sigma', sigma, 'Q', Q, 'Sig', Sig, 'QSigP.T', np.matmul(np.matmul(Q, Sig), P))
```

## Линейная регрессия
В некотором эксперименте измерялись значения величин $g_1$, $g_2$, $g_3$ и $H$:
$$
\begin{matrix}
g_1 & 0.12 & 0.15 & 0.9 & 0.8\\
g_2 & 2.4  & 1.8  & 3.2 & 3.6\\
g_3 & 1.1  & 1.2  & 1.3 & 1.4\\
H   & 5.1  & 6.2  & 5.5 & 4.1
\end{matrix}
$$

Найти коэффициенты $a$, $b$, $c$ линейной регрессии $H = ag_1 + bg_2 +cg_3$.

Составим матрицу $A$ столбцы которой образуют значения $g_1$, $g_2$, $g_3$.
Также составим матрицу-столбец $H$ из значений $H$,
тогда
$$
\left[\begin{matrix}a\\b\\c\end{matrix}\right] = A^+H
$$


```python
A = np.array([[ 0.12, 0.15, 0.9, 0.8], 
               [ 2.4, 1.8, 3.2, 3.6],
               [1.1, 1.2, 1.3, 1.4]]).transpose()
H1 = np.array([[5.1], [6.2], [5.5],  [4.1]])
res = np.matmul(np.linalg.pinv(A), H1)
a, b, c = res[:, 0]
display('A', A, 'H', H1, 'a', a, 'b', b, 'c', c)
```

Вычислим относительные отклонения экспериментальных данных от функции $H = ag_1 + bg_2 + cg_3$


```python
def Hfunc(g1, g2, g3):
    return a*g1 + b*g2 + c*g3
print(*[abs((Hfunc(*g) - H1[i][0])/H1[i][0]) for i, g in enumerate(A)])
```

60541b9c7393255e8754bd01e8a8d52469397cd7
ipynb
Jupyter Notebook
Applied_Alg_sem_4.ipynb
GalinaZh/Appl_alg2021
09761b56eb2bdfee4cd5f12cd96562ca146fcb15
MIT
Applied_Alg_sem_4.ipynb
GalinaZh/Appl_alg2021
09761b56eb2bdfee4cd5f12cd96562ca146fcb15
MIT
Applied_Alg_sem_4.ipynb
GalinaZh/Appl_alg2021
09761b56eb2bdfee4cd5f12cd96562ca146fcb15
MIT
Qwen/Qwen-72B
1. YES
2. YES