# Занятие 2
# Прикладная алгебра и численные методы
## Псевдообратная матрица, скелетное разложение
numpy:
https://numpy.org/doc/stable/reference/routines.linalg.html

scipy:
https://docs.scipy.org/doc/scipy/reference/linalg.html


```python
import numpy as np
import scipy.linalg
import sympy
```

## Псевдообратная матрица
В numpy псевдообратная матрица находится с помощью np.linalg.pinv
## Пример 1
Найти псевдообратную матрицу к матрице
$$
\left(
\begin{matrix}
1 & 2 & 3\\
4 & 5 & 6
\end{matrix}
\right)
$$


```python
A = np.arange(1, 7).reshape((2, 3))
A_pinv = np.linalg.pinv(A)
display(A, A_pinv)
np.allclose(A, np.dot(A, np.dot(A_pinv, A)))
np.allclose(A_pinv, np.dot(A_pinv, np.dot(A, A_pinv)))
```


    array([[1, 2, 3],
           [4, 5, 6]])



    array([[-0.94444444,  0.44444444],
           [-0.11111111,  0.11111111],
           [ 0.72222222, -0.22222222]])





    True



## В scipy тоже есть модуль linalg, 
в нем есть pinv:


```python
A = np.arange(1, 7).reshape((2, 3))
A_pinv = scipy.linalg.pinv(A)
display(A, A_pinv)
np.allclose(A, np.dot(A, np.dot(A_pinv, A)))
np.allclose(A_pinv, np.dot(A_pinv, np.dot(A, A_pinv)))
```


    array([[1, 2, 3],
           [4, 5, 6]])



    array([[-0.94444444,  0.44444444],
           [-0.11111111,  0.11111111],
           [ 0.72222222, -0.22222222]])





    True



## В пакете символьной математики sympy 
тоже есть псевдообратная матрица:


```python
As = sympy.Matrix(A)
As_pinv = As.pinv()
display(As_pinv, As_pinv.evalf(8))
```


$\displaystyle \left[\begin{matrix}- \frac{17}{18} & \frac{4}{9}\\- \frac{1}{9} & \frac{1}{9}\\\frac{13}{18} & - \frac{2}{9}\end{matrix}\right]$



$\displaystyle \left[\begin{matrix}-0.94444444 & 0.44444444\\-0.11111111 & 0.11111111\\0.72222222 & -0.22222222\end{matrix}\right]$



```python
np.array(As)
```




    array([[1, 2, 3],
           [4, 5, 6]], dtype=object)



## Скелетное разложение

Построим скелетное разложение для матрицы 
$$
\left(
\begin{matrix}
2 & 3 & -1 & 4\\
-1 & 4 & 1 & 4\\
1 & 7 & 0 & 8
\end{matrix}
\right)
$$
$B$ матрица полного столбцового ранга, $C$-строчного,
тогда псевдообратная матрица находится по формуле
$$
A^+ = C^T(C\cdot C^T)^{-1}(B^T\cdot B)^{-1}B^T
$$


```python
A = sympy.Matrix([[2, 3, -1, 4], [-1, 4, 1, 4], [1, 7, 0, 8]])
A_rref = A.rref()
display(*A_rref)
```


$\displaystyle \left[\begin{matrix}1 & 0 & - \frac{7}{11} & \frac{4}{11}\\0 & 1 & \frac{1}{11} & \frac{12}{11}\\0 & 0 & 0 & 0\end{matrix}\right]$



    (0, 1)



```python
cols = A_rref[1]
k = len(cols)
B = A[:, cols]
B
```




$\displaystyle \left[\begin{matrix}2 & 3\\-1 & 4\\1 & 7\end{matrix}\right]$




```python
C = A_rref[0][:k, :]
C
```




$\displaystyle \left[\begin{matrix}1 & 0 & - \frac{7}{11} & \frac{4}{11}\\0 & 1 & \frac{1}{11} & \frac{12}{11}\end{matrix}\right]$




```python
A_pinv_my = C.T*(C*C.T)**(-1)*(B.T*B)**(-1)*B.T
display(A, B, C, A_pinv_my, A.pinv() == A_pinv_my)
```


$\displaystyle \left[\begin{matrix}2 & 3 & -1 & 4\\-1 & 4 & 1 & 4\\1 & 7 & 0 & 8\end{matrix}\right]$



$\displaystyle \left[\begin{matrix}2 & 3\\-1 & 4\\1 & 7\end{matrix}\right]$



$\displaystyle \left[\begin{matrix}1 & 0 & - \frac{7}{11} & \frac{4}{11}\\0 & 1 & \frac{1}{11} & \frac{12}{11}\end{matrix}\right]$



$\displaystyle \left[\begin{matrix}\frac{266}{1185} & - \frac{253}{1185} & \frac{13}{1185}\\- \frac{41}{1185} & \frac{88}{1185} & \frac{47}{1185}\\- \frac{173}{1185} & \frac{169}{1185} & - \frac{4}{1185}\\\frac{52}{1185} & \frac{4}{1185} & \frac{56}{1185}\end{matrix}\right]$



    True


## Функции для работы с array.

https://numpy.org/doc/stable/reference/routines.array-manipulation.html

### Копирование
При копировании с помощью присваивания получается новый указатель на тот же объект, а не физически независимая копия. Но при использовании операции умножения на число результат другой, сравните:


```python
ar1 = np.array([[1, 2, 3, 4], [5, 6, 7, 8]])
ar2 = 2*ar1
display(ar1, ar2)
ar1[0, 0] = -1
ar2[1, 1] = 100
display(ar1, ar2)
```


    array([[1, 2, 3, 4],
           [5, 6, 7, 8]])



    array([[ 2,  4,  6,  8],
           [10, 12, 14, 16]])



    array([[-1,  2,  3,  4],
           [ 5,  6,  7,  8]])



    array([[  2,   4,   6,   8],
           [ 10, 100,  14,  16]])


ar1 и ar2 изменяются независимо

Теперь не будем умножать:


```python
ar1 = np.array([[1, 2, 3, 4], [5, 6, 7, 8]])
ar2 = ar1
display(ar1, ar2)
ar1[0, 0] = -1
ar2[1, 1] = 100
display(ar1, ar2)
```


    array([[1, 2, 3, 4],
           [5, 6, 7, 8]])



    array([[1, 2, 3, 4],
           [5, 6, 7, 8]])



    array([[ -1,   2,   3,   4],
           [  5, 100,   7,   8]])



    array([[ -1,   2,   3,   4],
           [  5, 100,   7,   8]])


Видно, что ar2 - это новый указатель на ar1.

Умножим на 1:


```python
ar1 = np.array([[1, 2, 3, 4], [5, 6, 7, 8]])
ar2 = 1*ar1
display(ar1, ar2)
ar1[0, 0] = -1
ar2[1, 1] = 100
display(ar1, ar2)
```


    array([[1, 2, 3, 4],
           [5, 6, 7, 8]])



    array([[1, 2, 3, 4],
           [5, 6, 7, 8]])



    array([[-1,  2,  3,  4],
           [ 5,  6,  7,  8]])



    array([[  1,   2,   3,   4],
           [  5, 100,   7,   8]])


Этот способ очень не рекомендуется применять в групповом или развивающемся проекте,
поскольку велики шансы, что кто-нибудь решит убрать "ненужную" единицу и все сломается.

Попробуйте все то же самое проделать с списком list. Тут умножение не помогает:


```python
list1 = [[1, 2, 3, 4], [5, 6, 7, 8]]
list2 = 1*list1
display(list1, list2)
list1[0][0] = -1
list2[1][1] = 100
display(list1, list2)
```


    [[1, 2, 3, 4], [5, 6, 7, 8]]



    [[1, 2, 3, 4], [5, 6, 7, 8]]



    [[-1, 2, 3, 4], [5, 100, 7, 8]]



    [[-1, 2, 3, 4], [5, 100, 7, 8]]


### Копирование с numpy.copyto
numpy.copyto(dst, src, casting='same_kind', where=True)

dst - куда (копировать)

src - что (копировать)

dst и src - 

Скопируем ar1 в ar2


```python
ar1 = np.array([[-1, 2, -3, 4], [5, -6, 7, -8]])
np.copyto(ar2, ar1)
display(ar1, ar2)
ar1[0, 0] = 555
ar2[1, 1] = 777
display(ar1, ar2)
```


    array([[-1,  2, -3,  4],
           [ 5, -6,  7, -8]])



    array([[-1,  2, -3,  4],
           [ 5, -6,  7, -8]])



    array([[555,   2,  -3,   4],
           [  5,  -6,   7,  -8]])



    array([[ -1,   2,  -3,   4],
           [  5, 777,   7,  -8]])


### numpy.asarray(a, dtype=None, order=None) 
возвращает array, полученный на основе аргумента a.


```python
a = [1, 2, 3]
b = np.asarray(a)
a, b
```




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



### numpy.concatenate((a1, a2, ...), axis=0, out=None)
соединяет аргументы в один array


```python
np.concatenate((ar1[:, :-1], [a]))
```




    array([[555,   2,  -3],
           [  5,  -6,   7],
           [  1,   2,   3]])



### numpy.stack(arrays, axis=0, out=None)
соединяет аргументы в один array, все аргументы должны быть одинаковых размерностей:


```python
ar3 = np.array([6, 7, 8])
ar4 = np.array([16, 0, -8])
np.stack((ar3, ar4))
```




    array([[ 6,  7,  8],
           [16,  0, -8]])



### numpy.vstack(tup)


```python
ar3 = np.array([[6, 3], [7, 4], [8, 5]])
ar4 = np.array([[1, 16], [2, 0], [3, -8]])
np.vstack((ar3, ar4))
```




    array([[ 6,  3],
           [ 7,  4],
           [ 8,  5],
           [ 1, 16],
           [ 2,  0],
           [ 3, -8]])



Сравните с numpy.stack


```python
np.stack((ar3, ar4))
```




    array([[[ 6,  3],
            [ 7,  4],
            [ 8,  5]],
    
           [[ 1, 16],
            [ 2,  0],
            [ 3, -8]]])



### numpy.hstack(tup)


```python
np.hstack((ar3, ar4))
```




    array([[ 6,  3,  1, 16],
           [ 7,  4,  2,  0],
           [ 8,  5,  3, -8]])



78e6453f88a431884bf852e76341deb5c8281489
ipynb
Jupyter Notebook
Applied_Alg_sem_2.ipynb
GalinaZh/Appl_alg2021
09761b56eb2bdfee4cd5f12cd96562ca146fcb15
MIT
Applied_Alg_sem_2.ipynb
GalinaZh/Appl_alg2021
09761b56eb2bdfee4cd5f12cd96562ca146fcb15
MIT
Applied_Alg_sem_2.ipynb
GalinaZh/Appl_alg2021
09761b56eb2bdfee4cd5f12cd96562ca146fcb15
MIT
Qwen/Qwen-72B
1. YES
2. YES