# NUME

Instructor: Keny **ORDAZ**
    
# Numerical Linear Algebra

$\newcommand{\vb}[1]{\boldsymbol{#1}}\newcommand{\RR}[1]{\mathbb{R}^{#1}}
\newcommand{\vt}[1]{\vb{#1}^\top}$


```python
import numpy as np
import scipy
import sympy as sp
import matplotlib.pyplot as plt
```

## Almacenamiento


```python
import numpy as np  # recomendado
# from numpy import *   # no es recomendado

from numpy import pi, cos, sin

a = np.array([0, 1, 2], dtype=np.float32)
b = np.array([0., 1, 2])
a.dtype, b.dtype
```




    (dtype('float32'), dtype('float64'))




```python
print(id(a))
```

    4582173456



```python
a +=  2*b
a
```




    array([0., 3., 6.], dtype=float32)




```python
print(id(a))
```

    4582173456



```python
a = a + 2 * b
```


```python
print(id(a))
```

    4800541104



```python
a, a.dtype
```




    (array([ 0.,  5., 10.]), dtype('float64'))




```python
z = np.zeros((5,1))
z, z.shape, z.dtype, len(z.shape), z[0,-1]
```




    (array([[0.],
            [0.],
            [0.],
            [0.],
            [0.]]),
     (5, 1),
     dtype('float64'),
     2,
     0.0)




```python
z = np.zeros((1, 5))
z, z.shape, z.dtype, len(z.shape), z[0,0]
```




    (array([[0., 0., 0., 0., 0.]]), (1, 5), dtype('float64'), 2, 0.0)




```python
z = np.zeros(5)
z, z.shape, z.dtype, len(z.shape), z[0]
```




    (array([0., 0., 0., 0., 0.]), (5,), dtype('float64'), 1, 0.0)




```python
nb = 100001
x = np.linspace(0, 2 * np.pi, nb, dtype=np.float32) # [0, 2π]
y = np.sin(x)
%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(x, y)
```


```python
x2 = np.zeros(nb)
x2 += x
x2.dtype, x.dtype
np.save('vec1001_32', x)
np.save('vec1001_64', x2)
```


```python
!ls -l *.npy
```

    -rw-r--r--@ 1 kordaz  staff     168 May 23 22:01 data-numpy-01.npy
    -rw-r--r--@ 1 kordaz  staff     216 May 23 22:23 data-numpy-02.npy
    -rw-r--r--@ 1 kordaz  staff     164 May 23 23:19 data-numpy-03.npy
    -rw-r--r--@ 1 kordaz  staff  400132 May 26 04:40 vec1001_32.npy
    -rw-r--r--@ 1 kordaz  staff  800136 May 26 04:40 vec1001_64.npy



```python
!dir
```

    zsh:1: command not found: dir



```python
xdisco = np.load('vec1001_64.npy')
```


```python
xdisco.dtype, xdisco.shape
```




    (dtype('float64'), (100001,))



## Representación de operaciones básicas

Sea $\alpha \in \mathbb{R},\quad \{\vb{A}, \vb{B}, \vb{C}\} \in \mathbb{R}^{m\times n},\quad \{\vb{D}\} \in \mathbb{R}^{m\times p},\quad \{\vb{E}\} \in \mathbb{R}^{p\times n}$

Transpuesta $\vt{M}$


```python
A = [[1, 2, 8], [3, 4, 5], [1, -1, 1]]
C = [[0, 0, 0], [0, 0, 0], [0, 0, 0]]

# transpuesta
for i in range(len(C)):
    for j in range(len(C[0])):
        C[i][j] = A[j][i]
```


```python
C
```




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




```python
import numpy as np
```


```python
matA = np.array(A)
matC = matA.T
print(matA)
matC
```

    [[ 1  2  8]
     [ 3  4  5]
     [ 1 -1  1]]





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




```python
B = np.arange(3*3).reshape((-1, 3))
B
C = matA + B
C
C = matA * B
C
```




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



$\vb{C}_{ij} = \vb{A}_{ij} \vb{B}_{ij}$


```python
A = matA
C = 0.4 * A
C
```




    array([[ 0.4,  0.8,  3.2],
           [ 1.2,  1.6,  2. ],
           [ 0.4, -0.4,  0.4]])




```python
m, p, n = 4, 5, 2
D = np.arange(m * p).reshape((m, p))
E = np.arange(p * n).reshape((p, n))
D, E
D @ E # O(m p n)
# NO factible: D * E
```




    array([[ 60,  70],
           [160, 195],
           [260, 320],
           [360, 445]])



Sea $\alpha, c \in \mathbb{R},\quad \{\vb{x}, \vb{y}, \vb{z}\} \in \mathbb{R}^{n}$


```python

```

#### Multiplicación 'vector' (elemento a elemento)


```python
x = np.arange(4, dtype=float); y = np.arange(4, 4*2, dtype=float); a = 1.5
x, y, a
```




    (array([0., 1., 2., 3.]), array([4., 5., 6., 7.]), 1.5)




```python
z = x + y
z
```




    array([ 4.,  6.,  8., 10.])



#### Producto punto

$\vb{x} \in \RR{m}, \vb{y} \in \RR{n}, \vb{x} \cdot \vb{y} = \vt{x} \vb{y} \in \RR{}$


```python
z = x @ y  # O(n)
(z, x.shape, y.shape)
#np.dot(x, y)
#x.dot(y)
x.reshape((-1,1)).T @ y.reshape((-1, 1)), x @ y
```




    (array([[38.]]), 38.0)




```python
z = x * y
z
```




    array([ 0.,  5., 12., 21.])



#### Actualización saxpy ("scalar a x plus y")


```python
y += x * a
y
```




    array([ 4. ,  6.5,  9. , 11.5])



#### Producto exterior

$\vb{x} \in \RR{m}, \vb{y} \in \RR{n}, \vb{x} \vt{y} \in \RR{m \times n}$


```python
x = np.array([1, 2, 3])
y = np.array([4, 5])
np.outer(x, y), x.reshape((-1,1)) @ y.reshape((-1, 1)).T
```




    (array([[ 4,  5],
            [ 8, 10],
            [12, 15]]),
     array([[ 4,  5],
            [ 8, 10],
            [12, 15]]))




```python
%timeit np.outer(x, y)
```

    3.16 µs ± 605 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)



```python
%timeit x.reshape((-1,1)) @ y.reshape((-1, 1)).T
```

    1.76 µs ± 224 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)


## Orden de complejidad


```python
x = np.linspace(0, 10000000000, 11)
y = 2*x**(2)
m = 0.1
y2 = m*x + 100000
y3 = 0.00000001 * x**3 - 456**3
plt.plot(x)
plt.plot(y, label='$\mathcal{O}(n^{2})$')
plt.plot(y2, label='$\mathcal{O}(n)$')
plt.plot(y3, label='$\mathcal{O}(n^3)$')
plt.legend()
```

Crecimiento con respecto al crecimiento de $n$


```python
x = np.arange(4, dtype=float); y = np.arange(4, 4*2, dtype=float); a = 1.5
x, y, a
```




    (array([0., 1., 2., 3.]), array([4., 5., 6., 7.]), 1.5)



### Producto punto


```python
c = 0
for i in range(x.shape[0]):
    c += x[i] * y[i]
c, x.shape
```




    (38.0, (4,))



c = c + x[i] * y[i]
(1 asig, 1 suma, 1 prod) * n
nb_ops_arit = 2 * n
O(n)


```python
A = np.arange(30).reshape((10, -1))
A
```




    array([[ 0,  1,  2],
           [ 3,  4,  5],
           [ 6,  7,  8],
           [ 9, 10, 11],
           [12, 13, 14],
           [15, 16, 17],
           [18, 19, 20],
           [21, 22, 23],
           [24, 25, 26],
           [27, 28, 29]])




```python
A[2]
```




    array([6, 7, 8])




```python
A[2,:]
```




    array([6, 7, 8])




```python
A[:, -1]
```




    array([ 2,  5,  8, 11, 14, 17, 20, 23, 26, 29])



ij -> ji... ijk jik, kji...

## Strassen Matrix Multiplication


```python
import time
import numpy as np
import numpy

def matrix_multiplication(A, B):
    C = np.zeros_like(A)
    m, _ = A.shape
    for i in range(m):
        for j in range(m):
            for k in range(m):
                C[i, j] += A[i, k] * B[k, j] 
    return C

def matrix_multiplication_dotprod(A, B):
    C = np.zeros_like(A)
    m, _ = A.shape
    for i in range(m):
        for j in range(m):
            C[i, j] += A[i, :] @ B[:, j] 
    return C

def strassen(A,B):
    n = A.shape[0]
    C = np.zeros((n, n))
    if n == 1:
        C[0][0] = A[0][0] * B[0][0]

    else:
        k = int(n/2) 

        A11,A21,A12,A22 = A[:k,:k], A[k:, :k], A[:k, k:], A[k:, k:]
        B11,B21,B12,B22 = B[:k,:k], B[k:, :k], B[:k, k:], B[k:, k:]

        S1 = B12 - B22
        S2 = A11 + A12
        S3 = A21 + A22
        S4 = B21 - B11
        S5 = A11 + A22
        S6 = B11 + B22
        S7 = A12 - A22
        S8 = B21 + B22
        S9 = A11 - A21
        S10= B11 + B12

        P1 = strassen(A11, S1)
        P2 = strassen(S2, B22)
        P3 = strassen(S3, B11)
        P4 = strassen(A22, S4)
        P5 = strassen(S5, S6)
        P6 = strassen(S7, S8)
        P7 = strassen(S9, S10)

        C[:k,:k] = P5 + P4 - P2 + P6
        C[:k, k:] = P1 + P2
        C[k:, :k] = P3 + P4
        C[k:, k:] = P5 + P1 - P3 - P7

    return C

"""
In linear algebra, the Strassen algorithm, named after Volker Strassen,
is an algorithm for matrix multiplication. It is faster than the standard matrix multiplication algorithm
and is useful in practice for large matrices, but would be slower than the fastest known algorithms
for extremely large matrices.
"""

if __name__ == '__main__':
    # Simple implementation needs n = 2**k
    k = 8
    n = 2 ** k
    A = np.random.rand(n,n) * 5
    B = np.random.rand(n,n) * 2 - 4

    tic = time.perf_counter_ns()
    D = A @ B
    toc = time.perf_counter_ns()
    print(f"Native matrix mult {' '*17} in {toc - tic:10d} nseconds")

    functions = [strassen, matrix_multiplication, matrix_multiplication_dotprod]

    for f in functions:
        tic = time.perf_counter_ns()
        C = f(A, B)
        toc = time.perf_counter_ns()
        print(f"Using {f.__name__:30} in {toc - tic:10d} nseconds, sum_error abs: ",
                np.sum(np.abs(C-D)))

    A = np.random.randint(20, size=(n,n))
    B = np.random.randint(20, size=(n,n))
    C = strassen(A, B)
    D = A @ B
    print(np.sum(np.abs(C-D)))

```

    Native matrix mult                   in    1350447 nseconds
    Using strassen                       in 22817056669 nseconds, sum_error abs:  5.299432359606726e-07
    Using matrix_multiplication          in 8723724799 nseconds, sum_error abs:  6.7700511863222346e-09
    Using matrix_multiplication_dotprod  in  122885372 nseconds, sum_error abs:  4.3129603000124916e-08
    0.0


# Resolución de sistemas lineales con métodos directos


ch03-c1


```python
import sympy as sp
import numpy as np

"""
m ecuaciones lineales y n unknowns
"""
def demo01(): # simbólico
    """

    Output
    ------
    {x0: 4, x1: -2, x2: 3}
    {x0: -1, x1: 2, x2: 0}
    """
    x = sp.symbols('x:3')  # x_0, x_1, x_2
    eq0 = 2 * x[0] + 4 * x[1] + 6 * x[2] - 18
    eq1 = 4 * x[0] + 5 * x[1] + 6 * x[2] - 24
    eq2 = 3 * x[0] + x[1] - 2 * x[2] - 4
    print(sp.solve([eq0, eq1, eq2], x))

    eq0 = 1 * x[0] +  2 * x[1] +  3 * x[2] -3
    eq1 = 4 * x[0] +  5 * x[1] +  6 * x[2] -6
    eq2 = 7 * x[0] +  8 * x[1] +  10 * x[2] -9
    print(sp.solve([eq0, eq1, eq2], x))
    print(type(x[0]))

def demo01symbolic():
    """

    Output
    ------
    FiniteSet((-1, 2, 0))
    """
    x, y, z = sp.symbols("x, y, z")
    A = sp.Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 10]])
    b = sp.Matrix([3, 6, 9])
    print(sp.linsolve((A, b), [x, y, z]))

def demo01numeric():
    """
    """
    x = np.array([-1, 2, 0])
    A = np.array([[1., 2, 3], [4, 5, 6], [7, 8, 10]])
    b = np.array([3., 6, 9])
    print("Sistema Ax =b, verificación: ", A @ x - b)

    print("\nOperaciones elementales")
    print(A)
    C = A.copy()
    C[[0, 1]] = C[[1, 0]]
    print(C, " intercambio de renglones")
    k = 2.05  # k != 0
    C[2, :] = k * C[2, :]
    print(C, " escalamiento de renglón")
    i, j, k = 0, 1, 3.0
    C[j, :] += k * C[i, :]
    print(C, "combinación lineal con un múltiplo de otro renglón")
    E = np.c_[A, b]
    print(f'E = \n{E}')
    i, j, k = 0, 1, -4.0
    E[j, :] += k * E[i, :]
    print(E)
    i, j, k = 0, 2, -7.0
    E[j, :] += k * E[i, :]
    print(E)
    i, j, k = 1, 2, -2.0
    E[j, :] += k * E[i, :]
    print(E)
    
    
    # Solución de sistemas mxm de ecuaciones lineales
    x_sol = np.linalg.solve(A, b)
    print("Solución numpy:", x_sol, "Solución esperada:", x, sep='\n')
    print(np.allclose(A @ x_sol, b))


if __name__ == '__main__':
    demo01symbolic()
    demo01()
    demo01numeric()


```

    FiniteSet((-1, 2, 0))
    {x0: 4, x1: -2, x2: 3}
    {x0: -1, x1: 2, x2: 0}
    <class 'sympy.core.symbol.Symbol'>
    Sistema Ax =b, verificación:  [0. 0. 0.]
    
    Operaciones elementales
    [[ 1.  2.  3.]
     [ 4.  5.  6.]
     [ 7.  8. 10.]]
    [[ 4.  5.  6.]
     [ 1.  2.  3.]
     [ 7.  8. 10.]]  intercambio de renglones
    [[ 4.    5.    6.  ]
     [ 1.    2.    3.  ]
     [14.35 16.4  20.5 ]]  escalamiento de renglón
    [[ 4.    5.    6.  ]
     [13.   17.   21.  ]
     [14.35 16.4  20.5 ]] combinación lineal con un múltiplo de otro renglón
    E = 
    [[ 1.  2.  3.  3.]
     [ 4.  5.  6.  6.]
     [ 7.  8. 10.  9.]]
    [[ 1.  2.  3.  3.]
     [ 0. -3. -6. -6.]
     [ 7.  8. 10.  9.]]
    [[  1.   2.   3.   3.]
     [  0.  -3.  -6.  -6.]
     [  0.  -6. -11. -12.]]
    [[ 1.  2.  3.  3.]
     [ 0. -3. -6. -6.]
     [ 0.  0.  1.  0.]]
    Solución numpy:
    [-1.00000000e+00  2.00000000e+00 -2.22044605e-16]
    Solución esperada:
    [-1  2  0]
    True



```python
### Numpy Bloques
```


```python
np.block([[np.eye(3)*5, np.zeros((3,2))],
         [np.zeros((2,3)), -2.5*np.eye(2)]])
```




    array([[ 5. ,  0. ,  0. ,  0. ,  0. ],
           [ 0. ,  5. ,  0. ,  0. ,  0. ],
           [ 0. ,  0. ,  5. ,  0. ,  0. ],
           [ 0. ,  0. ,  0. , -2.5, -0. ],
           [ 0. ,  0. ,  0. , -0. , -2.5]])



## Sistemas triangulares

### Gauss Elimination

The differences in the ranges of k in python and MATLAB reflect the native offsets used for arrays. In Python, all sequences have zero offset while in MATLAB the native offset is 1.


```python
function [x,det] = gaussElimin(a,b)
n = length(b);
for k = 1:n-1
    for i = k+1:n
        if a(i,k)~=0
            lam = a(i,k)/a(k,k);
            a(i,k+1:n) = a(i,k+1:n)-lam*a(k,K+1:n);
            b(i) = b(i)-lam*b(k);
        end
    end
end
det = prod(diag(a));
for k = n:-1:1
    b(k) = (b(k)-a(k,k+1:n)*b(k+1:n))/a(k,k);
end
x = b;
```



```python
from numpy import dot
from numpy import array

from book.gaussElimin import gaussElimin

a = np.array([[1.,1,1,1],[0,1,-1,-1],[0,0,-5,-8],[0,0,0,11]])

b = np.array([[10.],[-5],[-47],[44]])
print("Eliminación gaussiana\na y b son modificados al terminar el proceso")
print("Sistema 1\n", (x:=gaussElimin(a,b)), a@x - b)  

a = np.array([[1.,-2,2],[2,-3,-3],[4,1,6]])

b = np.array([[-2.],[4],[3]])

x_npy = np.linalg.solve(a, b)
print("Sistema 2\nSolución numpy:", x_npy, a @ x_npy - b, sep='\n')
print(np.allclose(np.dot(a, x_npy), b, 1e-10, 1e-10))
print(a)
print("Solución gaussElimin (llamada con copias):", (x:=gaussElimin(a.copy(), b.copy())), a@x - b, sep='\n')  
print(a)
```

    Eliminación gaussiana
    a y b son modificados al terminar el proceso
    Sistema 1
     [[1.]
     [2.]
     [3.]
     [4.]] [[  9.]
     [ -7.]
     [-50.]
     [ 40.]]
    Sistema 2
    Solución numpy:
    [[ 2.]
     [ 1.]
     [-1.]]
    [[0.]
     [0.]
     [0.]]
    True
    [[ 1. -2.  2.]
     [ 2. -3. -3.]
     [ 4.  1.  6.]]
    Solución gaussElimin (llamada con copias):
    [[ 2.]
     [ 1.]
     [-1.]]
    [[0.]
     [0.]
     [0.]]
    [[ 1. -2.  2.]
     [ 2. -3. -3.]
     [ 4.  1.  6.]]


## Factorización LU


```python
import numpy as np
from scipy.linalg import lu as LU
A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
p, l, u = LU(A)
res = np.allclose(A - p @ l @ u, np.zeros((4, 4)))

print("P L U == A?", res)
A = np.array([[ 1, 2, 3, 1],
    [ 2, 4, -1, 1],
    [ 1, 3, -1, 2],
    [ -3, -1, -3, 2]])
b = np.array([8, 1, 1, -7])
print("Matriz A", A, "vector b", b, sep='\n')

from scipy.linalg import lu_factor, lu_solve
lu, piv = lu_factor(A)
p, l, u = LU(A)
x = lu_solve((lu, piv), b)

x = lu_solve(lu_factor(A)b, b)


res = np.allclose(A@x, b)
with np.printoptions(precision=2, suppress=True, threshold=5):
    print("A x ≈ b?", res, "mat A", A, " P, L, U", p, l , u, sep='\n') 

```

    P L U == A? True
    Matriz A
    [[ 1  2  3  1]
     [ 2  4 -1  1]
     [ 1  3 -1  2]
     [-3 -1 -3  2]]
    vector b
    [ 8  1  1 -7]
    A x ≈ b?
    True
    mat A
    [[ 1  2  3  1]
     [ 2  4 -1  1]
     [ 1  3 -1  2]
     [-3 -1 -3  2]]
     P, L, U
    [[0. 0. 1. 0.]
     [0. 1. 0. 0.]
     [0. 0. 0. 1.]
     [1. 0. 0. 0.]]
    [[ 1.    0.    0.    0.  ]
     [-0.67  1.    0.    0.  ]
     [-0.33  0.5   1.    0.  ]
     [-0.33  0.8   0.11  1.  ]]
    [[-3.   -1.   -3.    2.  ]
     [ 0.    3.33 -3.    2.33]
     [ 0.    0.    3.5   0.5 ]
     [ 0.    0.    0.    0.74]]



```python
M = np.array([[3, 5], [6, 7]])
p, l, u = LU(M)
p, l, u

L = np.array([[1, 0], [2, 1]])
U = np.array([[3, 5], [0, -3]])
l, u, L, U, L @ U - M, p @ l @ u - M

#L + U - np.eye(2)
```




    (array([[1. , 0. ],
            [0.5, 1. ]]),
     array([[6. , 7. ],
            [0. , 1.5]]),
     array([[1, 0],
            [2, 1]]),
     array([[ 3,  5],
            [ 0, -3]]),
     array([[0, 0],
            [0, 0]]),
     array([[0., 0.],
            [0., 0.]]))




```python
import sympy as sp
A = sp.Matrix([[1, 2, 3],
               [4, 5, 6],
               [10, 11, 9]])
A
```




$\displaystyle \left[\begin{matrix}1 & 2 & 3\\4 & 5 & 6\\10 & 11 & 9\end{matrix}\right]$




```python
L, U, _ = A.LUdecomposition()
L, U, L @ U - A
```




    (Matrix([
     [ 1, 0, 0],
     [ 4, 1, 0],
     [10, 3, 1]]),
     Matrix([
     [1,  2,  3],
     [0, -3, -6],
     [0,  0, -3]]),
     Matrix([
     [0, 0, 0],
     [0, 0, 0],
     [0, 0, 0]]))




```python
M = sp.Matrix([[3, 5], [6, 7]])
lM, uM, _ = M.LUdecomposition()
lM, uM, lM @ uM - M
```




    (Matrix([
     [1, 0],
     [2, 1]]),
     Matrix([
     [3,  5],
     [0, -3]]),
     Matrix([
     [0, 0],
     [0, 0]]))



## Factorización LDMT

$$\vb{A} = \vb{L} \vb{D} \vt{M},  \vb{A} = \vb{L} \vb{U}$$

$$\vb{D} = [u_{ii}]$$

$$ \vt{M} = \vb{D}^{-1} \vb{U}$$




```python
import numpy as np
from scipy.linalg import lu as LU
A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
p, l, u = LU(A)
res = np.allclose(A - p @ l @ u, np.zeros((4, 4)))

d = np.diag(u)
mt = np.diag(d**(-1)) @ u

with np.printoptions(precision=2, suppress=True, threshold=5):
    print("P L U == A?", res, "mat A", A, " P, L, U", p, l , u, sep='\n') 
    print(d, mt,  p @ l @ np.diag(d) @ mt - A, sep='\n')
    
```

    P L U == A?
    True
    mat A
    [[2 5 8 7]
     [5 2 2 8]
     [7 5 6 6]
     [5 4 4 8]]
     P, L, U
    [[0. 1. 0. 0.]
     [0. 0. 0. 1.]
     [1. 0. 0. 0.]
     [0. 0. 1. 0.]]
    [[ 1.    0.    0.    0.  ]
     [ 0.29  1.    0.    0.  ]
     [ 0.71  0.12  1.    0.  ]
     [ 0.71 -0.44 -0.46  1.  ]]
    [[ 7.    5.    6.    6.  ]
     [ 0.    3.57  6.29  5.29]
     [ 0.    0.   -1.04  3.08]
     [ 0.    0.    0.    7.46]]
    [ 7.    3.57 -1.04  7.46]
    [[ 1.    0.71  0.86  0.86]
     [ 0.    1.    1.76  1.48]
     [ 0.    0.    1.   -2.96]
     [ 0.    0.    0.    1.  ]]
    [[ 0.  0. -0. -0.]
     [-0. -0.  0.  0.]
     [ 0. -0.  0.  0.]
     [-0. -0.  0.  0.]]


## Factorización LDLT

$\vb{A} = \vb{L} \vb{D} \vt{L}, \vb{A} = \vt{A}$


```python
A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
A = A + A.T # A es simétrica ahora
```


```python
lu, d, perm = scipy.linalg.ldl(A, lower=0)
lu, lu[perm], d, lu @ d @ lu.T - A
```




    (array([[ 1.        , -0.32758621,  1.30434783,  0.75      ],
            [ 0.        ,  1.        , -0.08695652,  0.75      ],
            [ 0.        ,  0.        ,  1.        ,  0.625     ],
            [ 0.        ,  0.        ,  0.        ,  1.        ]]),
     array([[ 1.        , -0.32758621,  1.30434783,  0.75      ],
            [ 0.        ,  1.        , -0.08695652,  0.75      ],
            [ 0.        ,  0.        ,  1.        ,  0.625     ],
            [ 0.        ,  0.        ,  0.        ,  1.        ]]),
     array([[-14.24137931,   0.        ,   0.        ,   0.        ],
            [  0.        ,  -5.04347826,   0.        ,   0.        ],
            [  0.        ,   0.        ,   5.75      ,   0.        ],
            [  0.        ,   0.        ,   0.        ,  16.        ]]),
     array([[0., 0., 0., 0.],
            [0., 0., 0., 0.],
            [0., 0., 0., 0.],
            [0., 0., 0., 0.]]))




```python
lu, d, perm = scipy.linalg.ldl(A, lower=True)
lu, lu[perm], d, lu @ d @ lu.T - A
```




    (array([[ 1.25      , -0.06521739,  1.        ,  0.        ],
            [ 0.58333333,  0.80434783, -0.11176471,  1.        ],
            [ 1.        ,  0.        ,  0.        ,  0.        ],
            [ 0.83333333,  1.        ,  0.        ,  0.        ]]),
     array([[ 1.        ,  0.        ,  0.        ,  0.        ],
            [ 0.83333333,  1.        ,  0.        ,  0.        ],
            [ 1.25      , -0.06521739,  1.        ,  0.        ],
            [ 0.58333333,  0.80434783, -0.11176471,  1.        ]]),
     array([[ 12.        ,   0.        ,   0.        ,   0.        ],
            [  0.        ,   7.66666667,   0.        ,   0.        ],
            [  0.        ,   0.        , -14.7826087 ,   0.        ],
            [  0.        ,   0.        ,   0.        ,  -4.85882353]]),
     array([[ 1.77635684e-15,  0.00000000e+00,  0.00000000e+00,
             -1.77635684e-15],
            [ 0.00000000e+00, -8.88178420e-16, -8.88178420e-16,
              0.00000000e+00],
            [ 0.00000000e+00, -8.88178420e-16,  0.00000000e+00,
              0.00000000e+00],
            [ 0.00000000e+00, -1.77635684e-15,  0.00000000e+00,
              0.00000000e+00]]))




```python
A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
Asim = 0.5 * (A + A.T) # A es simétrica ahora
Aanti = 0.5 * (A - A.T)
Asim, Aanti
Asim == Asim.T
Aanti == -Aanti.T
A == Asim + Aanti
```




    array([[ True,  True,  True,  True],
           [ True,  True,  True,  True],
           [ True,  True,  True,  True],
           [ True,  True,  True,  True]])



## Factorización de Cholesky

Simétrica, positiva definida

$$\vb{A} \in \RR{n \times n}, \vb{A} = \vt{A}, \vt{x}\vb{A}\vb{x} > 0 \quad \forall \vb{x} \neq \vb{0} $$



$$\vt{x}\vb{x} > 0\quad \forall \vb{x} \neq \vb{0}$$

Positiva semidefinida

$$\vt{x}\vb{A}\vb{x} \ge 0 \quad \forall \vb{x} \neq \vb{0}$$


```python
A = np.array([[2., 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
A = A @ A.T
A
```




    array([[142.,  92., 129., 118.],
           [ 92.,  97., 105., 105.],
           [129., 105., 146., 127.],
           [118., 105., 127., 121.]])




```python
w, v = np.linalg.eig(A)
with np.printoptions(precision=2, suppress=True, threshold=5):
    print(f'w: {w}', v, sep='\n')
```

    w: [468.08  27.22  10.42   0.28]
    [[ 0.52  0.74 -0.39 -0.18]
     [ 0.43 -0.63 -0.38 -0.53]
     [ 0.54  0.01  0.82 -0.17]
     [ 0.5  -0.24 -0.17  0.81]]


$$\vb{A} \vb{x} = \lambda \vb{x}$$


```python
with np.printoptions(precision=3, suppress=True, threshold=5):
    for i in range(len(w)):
        print(np.allclose(A @ v[:,i],  w[i] * v[:,i]))
```

    True
    True
    True
    True



```python
L = np.linalg.cholesky(A)
L, L @ L.T - A, L.T @ L - A
```

    True





    (array([[11.91637529,  0.        ,  0.        ,  0.        ],
            [ 7.7204685 ,  6.11509331,  0.        ,  0.        ],
            [10.82543952,  3.50322295,  4.06660646,  0.        ],
            [ 9.90234003,  4.66866069,  0.84774803,  0.65466982]]),
     array([[ 0.00000000e+00,  0.00000000e+00, -2.84217094e-14,
              0.00000000e+00],
            [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
              0.00000000e+00],
            [-2.84217094e-14,  0.00000000e+00,  0.00000000e+00,
              0.00000000e+00],
            [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
              0.00000000e+00]]),
     array([[ 274.85211268,   39.36597907,  -76.58250839, -111.5172368 ],
            [  39.36597907,  -25.53667011,  -86.795923  , -101.94356873],
            [ -76.58250839,  -86.795923  , -128.74403514, -126.44500495],
            [-111.5172368 , -101.94356873, -126.44500495, -120.57140742]]))




```python
A = np.array([[3, -6], [4, -8], [0, 1]])
b = np.array([1, 7, 2])

C = A.T @ A
C, (D:=A @ A.T), np.all(np.linalg.eigvals(C) >0), \
np.all(np.linalg.eigvals(D)>0)

L = np.linalg.cholesky(C)
d = A.T @ b # d: [ 31 -60]  TODO: corregir lecture notes

z = np.zeros(2)
x = np.zeros(2)
z[0] = d[0] / L[0,0] 
z[1] = (d[1] - L[1, 0] * z[0]) / L[1, 1]
x[1] = z[1] / L[1, 1]
x[0] = (z[0] - L[1, 0] * x[1]) / L[0, 0]

print(f'L\n{L}\nd: {d}\nz: {z}\n    x: {x} ')

c, low = scipy.linalg.cho_factor(C)
x_py = scipy.linalg.cho_solve((c, low), d)
from scipy.optimize import lsq_linear
res = lsq_linear(A, b)
np.allclose(res.x, x), x, x_py
res
```

    L
    [[  5.   0.]
     [-10.   1.]]
    d: [ 31 -60]
    z: [6.2 2. ]
        x: [5.24 2.  ] 





     active_mask: array([0., 0.])
            cost: 5.78
             fun: array([ 2.72, -2.04,  0.  ])
         message: 'The unconstrained solution is optimal.'
             nit: 0
      optimality: 7.815970093361102e-14
          status: 3
         success: True
               x: array([5.24, 2.  ])



#### Sensibilidad numérica


```python
%timeit np.linalg.cond(A.T @ A)
```

    30.7 µs ± 3.17 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)



```python
def cond(X):
    return np.linalg.norm(X) * np.linalg.norm(np.linalg.inv(X))
```


```python
%timeit cond(A.T @ A)
```

    22.6 µs ± 623 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)



```python
A = np.array([[1, 1], [1e-3, 0], [0, 1e-3]])
A, np.linalg.cond(A)
```




    (array([[1.   , 1.   ],
            [0.001, 0.   ],
            [0.   , 0.001]]),
     1414.2139159264418)




```python
U, s, V = scipy.linalg.svd(A)
S = np.zeros((3, 2))
S[:2, :2] = np.diag(s)
with np.printoptions(precision=3, suppress=True, threshold=5):
    print(U, s, V, U.shape, s.shape, V.shape, s[0] / s[1], U @ S @ V.T - A)
```

    [[-1.    -0.    -0.001]
     [-0.    -0.707  0.707]
     [-0.     0.707  0.707]] [1.414 0.001] [[-0.707 -0.707]
     [-0.707  0.707]] (3, 3) (2,) (2, 2) 1414.2139159264416 [[ 0.  0.]
     [-0. -0.]
     [-0.  0.]]



```python
# A x = b, A x = d
A = np.array([[1,1], [1,1.001]])
b = np.array([2, 2.])
d = np.array([2, 2.001])


np.linalg.solve(A, b), np.linalg.solve(A, d), np.linalg.cond(A)
```




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




```python
B = np.array([[0.001, 1],[1,1]])
e = np.array([1, 2.])
f = np.array([1-0.01, 2.001])
np.linalg.solve(B, e), np.linalg.solve(B, f), np.linalg.cond(B)
```




    (array([1.001001, 0.998999]),
     array([1.01201201, 0.98898799]),
     2.621550332216381)




```python

```

## Factorización QR

Factor the matrix a as qr, where q is orthonormal and r is upper-triangular.

$$ \vt{Q} = \vb{Q}^{-1} \leftarrow \vt{Q}\vb{Q} = \vb{I}$$


```python
a = np.random.randn(6, 3)
q, r = np.linalg.qr(a, 'complete')

with np.printoptions(precision=3, suppress=True, threshold=5):
    print(np.allclose(a, q @ r), q, r, sep='\n' ) # a does equal qr

```

    True
    [[-0.366 -0.119  0.387 -0.742 -0.364  0.139]
     [ 0.815 -0.138  0.291 -0.349  0.129 -0.305]
     [-0.007 -0.749 -0.13   0.25  -0.514 -0.309]
     [-0.135 -0.154  0.843  0.453  0.188  0.076]
     [-0.293 -0.513 -0.164 -0.244  0.742 -0.114]
     [ 0.311 -0.345 -0.099  0.013  0.002  0.88 ]]
    [[-3.061  0.264 -1.652]
     [ 0.    -2.697  1.647]
     [ 0.     0.     0.776]
     [ 0.     0.     0.   ]
     [ 0.     0.     0.   ]
     [ 0.     0.     0.   ]]



```python
np.allclose(q @ q.T, np.eye(6)), np.linalg.norm(q[0])
```




    (True, 0.9999999999999999)




```python
for qi in q:
    print(np.allclose(np.linalg.norm(qi), 1.0))
```

    True
    True
    True
    True
    True
    True



```python
t = np.pi / 5
R = np.array([[np.cos(t), -np.sin(t)], [np.sin(t), np.cos(t)]])
R @ R.T , np.linalg.det(R)
```




    (array([[1.00000000e+00, 1.26599665e-17],
            [1.26599665e-17, 1.00000000e+00]]),
     1.0000000000000002)




```python
v = np.array([15., -1., 0])

def build_householder_projector(v):
    #P = I - 2 * v v^T / v^T v
    return np.eye(len(v)) - 2 * np.outer(v, v) / np.dot(v, v)

P = build_householder_projector(v)
P, np.allclose(P @ P.T, np.eye(3)), np.allclose(np.abs(np.linalg.det(P)), 1)
```




    (array([[-0.99115044,  0.13274336,  0.        ],
            [ 0.13274336,  0.99115044,  0.        ],
            [ 0.        ,  0.        ,  1.        ]]),
     True,
     True)




```python
fig = plt.figure(figsize=(12, 8))
plt.gca().set_aspect('equal')
ax = fig.add_subplot(projection='3d')

p = np.array([0., 0, 0])

pp = np.array([p, p + v])
pq = np.array([p, [-v[1], v[0], 0]])

plt.plot(*pp[:,:2].T, 'k*:')
plt.plot(*pq[:,:2].T, 'r*:')

p0 = np.array([6, 1, 0.])
p1 = np.array([3, 8, 0.])

s = np.array([p0, p1])
plt.plot(*s[:,:2].T, 'bo-')

sp = P @ s.T
sp, P.shape, s.shape, s
plt.plot(*sp[:2,:], 'bo-')

q0 = np.array([4.5, 4.5, 0])

q0p = P @ q0
plt.plot(*q0[:2], 'g*-')
plt.plot(*q0p[:2], 'go-')
```


```python
v = np.array([15., -1., 1])
P = build_householder_projector(v)


fig = plt.figure(figsize=(12, 8))
plt.gca().set_aspect('equal')
ax = fig.add_subplot(projection='3d')

p = np.array([0., 0, 0])

pp = np.array([p, p + v])
pq = np.array([p, [-v[1], v[0], 0]])

plt.plot(*pp.T, 'k*:')
plt.plot(*pq.T, 'r*:')

p0 = np.array([6, 1, 0.])
p1 = np.array([3, 8, 0.])

s = np.array([p0, p1])
plt.plot(*s.T, 'bo-')

sp = P @ s.T
sp, P.shape, s.shape, s
plt.plot(*sp, 'bo-')

q0 = np.array([4.5, 4.5, 0])

q0p = P @ q0
plt.plot(*q0, 'g*-')
plt.plot(*q0p, 'go-')
```


```python
!ls book/*.py
```

    book/LUdecomp.py       book/example8_8.py     book/newtonRaphson2.py
    book/LUdecomp3.py      book/gaussElimin.py    book/plotPoly.py
    book/LUdecomp5.py      book/gaussNodes.py     book/polyFit.py
    book/LUpivot.py        book/gaussPivot.py     book/polyRoots.py
    book/bisection.py      book/gaussQuad.py      book/powell.py
    book/bulStoer.py       book/gaussQuad2.py     book/printSoln.py
    book/choleski.py       book/gaussSeidel.py    book/rational.py
    book/conjGrad.py       book/gerschgorin.py    book/ridder.py
    book/cubicSpline.py    book/goldSearch.py     book/romberg.py
    book/downhill.py       book/householder.py    book/rootsearch.py
    book/eigenvals3.py     book/inversePower.py   book/run_kut4.py
    book/error.py          book/inversePower3.py  book/run_kut5.py
    book/euler.py          book/inversePower5.py  book/sortJacobi.py
    book/evalPoly.py       book/jacobi.py         book/stdForm.py
    book/example8_1.py     book/lamRange.py       book/sturmSeq.py
    book/example8_3.py     book/linInterp.py      book/swap.py
    book/example8_4.py     book/midpoint.py       book/trapezoid.py
    book/example8_5.py     book/neville.py        book/triangleQuad.py
    book/example8_6.py     book/newtonPoly.py
    book/example8_7.py     book/newtonRaphson.py



```python
!cat book/householder.py
```

    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    


# Resolución de sistemas lineales con métodos iterativos


```python
import sympy as sp
from sympy.abc import a, b, c, d, u, w



```


```python
p0, pu0, pm, rhoH = (np.array([ 1.78201305, -0.4539905, 0 ]),
 np.array([0.907981  , 0.89100652, 0]),
 np.array([1.20002352, 0.24853301, 0]),
 np.array([0.64773555, 0.13415042, 0]))
```


```python
a, b, c, d = p0, pu0, pm, rhoH
```


```python
c @ (np.cross(a, d))
```




    0.0




```python
np.cross(a,b)
```




    array([-0.        ,  0.        ,  1.99999999])




```python
c
```




    array([1.20002352, 0.24853301, 0.        ])




```python
b @ (np.cross(a, d))
```




    0.0




```python
u = np.linalg.norm(c) / np.linalg.norm(b)
```


```python
p2 = a + u * b
```


```python
p2
```




    array([2.65670264, 0.40434698, 0.        ])




```python

```


```python

```


```python

```


```python
import numpy as np

def solve_sor(A, b, w, x):
    x = np.asarray(x)
    n = b.shape[0]
    xold = x.copy()
    convergence = False
    while not convergence:
        for i in range(n):
            s = 0.0
            for j in range(n):
                if j != i:
                    s += A[i,j]*x[j]
            suma = np.sum(A[i,:] @ x) - A[i,i]*x[i]
            #print(f'{suma=} \t\t {s=}\t{A[i,:] @ x}')
            x[i] = x[i] + w * ( (b[i] - s) / A[i, i] - x[i])
            # print(f'\t{x=}')
        # print(f'{x=}')
        convergence = np.allclose(A @ x, b, 1e-14) # or np.allclose(xold, x, 1e-14)
        xold = x.copy()
    return  x

def solve_jacobi_no(A, b, w, x):
    x = np.asarray(x)
    L, D, U = -np.tril(A, -1), np.diagonal(A), -np.triu(A, 1)
    assert (A == -L + np.diag(D) + -U).all()
    n = b.shape[0]
    xold = x.copy()
    convergence = False 
    k = 0
    while not convergence:
        k += 1
        x = ((L + U )@x - b)
        print(f'1 {x=}')
        x = x / D
        print(f'2 {x=}')
        convergence = k > 10 or np.allclose(A @ x, b, 1e-14) # or np.allclose(xold, x, 1e-14)
        xold = x.copy()
    return  x

def solve_jacobi(A, b, w, x, max_iter=1000):
    x = np.asarray(x)
    L, D, U = -np.tril(A, -1), np.diagonal(A), -np.triu(A, 1)
    assert (A == -L + np.diag(D) + -U).all()
    n = b.shape[0]
    xold = x.copy()
    convergence = False 
    k = 0
    while not convergence:
        k += 1
        for i in range(n):
            s = 0.0
            for j in range(n):
                if j != i:
                    s += A[i,j]*xold[j]
            suma = np.sum(A[i,:] @ x) - A[i,i]*x[i]
            #print(f'{suma=} \t\t {s=}\t{A[i,:] @ x}')
            x[i] = xold[i] + ( (b[i] - s) / A[i, i] - xold[i])
            # print(f'\t{x=}')
        # print(f'{x=}')
        convergence = k > max_iter or np.allclose(A @ x, b, 1e-14) # or np.allclose(xold, x, 1e-14)
        xold = x.copy()
    return  x


A = np.array([
    [4, -1, -6, 0.],
    [-5, -4, 10, 8],
    [0, 9, 4, -2],
    [1, 0, -7, 5]
    ])
b = np.array([2., 21, -12, -6]) 

E2_A = np.array([[10, -1, 2,0], [-1, 11, -1, 3], [2, -1, 10, -1], [3, -1, 8, 15]])
E2_b = np.array([6., 25, -11, 15])

for solver in [solve_sor, solve_jacobi]:
    x = solver(A, b, 0.5, [0.0, 0, 0, 0])
    print(x, A@x - b)

y = solve_jacobi(E2_A, E2_b, 0. , np.zeros_like(E2_b))
print(f'sol: {y}\n{b - A @ y}')

```


```python

```

# Matrices especiales


```python
import struct
import numpy as np
import scipy
from scipy.sparse import diags


def sliding_windows(a, W):
    """

    Parameters
    ----------

    See also
    --------

    Examples
    --------


    """
    a = np.asarray(a)
    p = np.zeros(W - 1, dtype=a.dtype)
    b = np.concatenate((p, a, p))
    s = b.strides[0]
    strided = np.lib.stride_tricks.as_strided
    return strided(b[W-1:], shape=(W, len(a)+W-1), strides=(-s, s))


def build_tridiag(a10, a00, a01, n=10):
    o = n - 1
    return np.diag([a00]*n) + np.diag([a10]*o, k=-1) + np.diag([a01]*o, k=+1)

def inverse_tridiag(n=10):
    temp = np.zeros((n, n))
    o = n + 1
    for j in range(n):
        for i in range(j + 1):
            #print(i <= j)
            temp[i, j] = temp[j, i] = (- i - 1) * (o - j - 1) / (o)
    return temp

if __name__ == "__main__":
    print("Construcción de matrices especiales")
    a = [1, -2., 1]
    A = sliding_windows(a, W=5)

    nb = 6
    B = diags(a, [0, 1, 2], shape=(nb, nb)).toarray()
    C = build_tridiag(*a, n=nb)
    D = inverse_tridiag(n=nb)

    for X in [A, B, C, D]:
        print(X)
    Iden2 =  C @ D
    
    print("Las matrices tridiagonales tienen simplificada la inversa y el determinante")
    for I in [Iden2]:
        print(np.allclose(I, np.eye(nb)))
        with np.printoptions(precision=3, suppress=True, threshold=5):
            print(I)
```

    Construcción de matrices especiales
    [[ 1. -2.  1.  0.  0.  0.  0.]
     [ 0.  1. -2.  1.  0.  0.  0.]
     [ 0.  0.  1. -2.  1.  0.  0.]
     [ 0.  0.  0.  1. -2.  1.  0.]
     [ 0.  0.  0.  0.  1. -2.  1.]]
    [[ 1. -2.  1.  0.  0.  0.]
     [ 0.  1. -2.  1.  0.  0.]
     [ 0.  0.  1. -2.  1.  0.]
     [ 0.  0.  0.  1. -2.  1.]
     [ 0.  0.  0.  0.  1. -2.]
     [ 0.  0.  0.  0.  0.  1.]]
    [[-2.  1.  0.  0.  0.  0.]
     [ 1. -2.  1.  0.  0.  0.]
     [ 0.  1. -2.  1.  0.  0.]
     [ 0.  0.  1. -2.  1.  0.]
     [ 0.  0.  0.  1. -2.  1.]
     [ 0.  0.  0.  0.  1. -2.]]
    [[-0.85714286 -0.71428571 -0.57142857 -0.42857143 -0.28571429 -0.14285714]
     [-0.71428571 -1.42857143 -1.14285714 -0.85714286 -0.57142857 -0.28571429]
     [-0.57142857 -1.14285714 -1.71428571 -1.28571429 -0.85714286 -0.42857143]
     [-0.42857143 -0.85714286 -1.28571429 -1.71428571 -1.14285714 -0.57142857]
     [-0.28571429 -0.57142857 -0.85714286 -1.14285714 -1.42857143 -0.71428571]
     [-0.14285714 -0.28571429 -0.42857143 -0.57142857 -0.71428571 -0.85714286]]
    Las matrices tridiagonales tienen simplificada la inversa y el determinante
    True
    [[ 1.  0.  0.  0.  0.  0.]
     [ 0.  1.  0. -0.  0.  0.]
     [-0. -0.  1.  0.  0.  0.]
     [ 0.  0.  0.  1. -0. -0.]
     [ 0.  0. -0.  0.  1.  0.]
     [ 0.  0.  0.  0.  0.  1.]]



```python
C
```


```python
np.linalg.det(C)
```


```python
p, l ,u = LU(C)
```


```python
p, l, u
```


```python
np.linalg.det(u)
```


```python
A = np.array([[1., 2, 3], [4, 5, 6], [7, 8, 10]])
b = np.array([3., 6, 9])

np.block([[A, b.reshape((-1, 1))]]), A.shape, b.shape
np.c_[A, b]
```




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




```python

```

5a8011f0c34c1214a630f7ecafad4271ec0f6c8f
ipynb
Jupyter Notebook
src/nume-s2-p1-float-numelinalg.ipynb
krontzo/nume.py
9d1e576fb3474333a8e2cf4f26f4236ee4f9deea
MIT
src/nume-s2-p1-float-numelinalg.ipynb
krontzo/nume.py
9d1e576fb3474333a8e2cf4f26f4236ee4f9deea
MIT
src/nume-s2-p1-float-numelinalg.ipynb
krontzo/nume.py
9d1e576fb3474333a8e2cf4f26f4236ee4f9deea
MIT
Qwen/Qwen-72B
1. YES
2. YES