# Mínimos Cuadrados

## Descomposición $QR$


```python
import numpy as np
import scipy.linalg as spla
import matplotlib.pyplot as plt
```


```python
def plotVector(A, Q):
    plt.figure(figsize=(6, 6))
    origin = [0], [0] # origin point
    plt.quiver(*origin, A[:,0], A[:,1], angles='xy', scale_units='xy', scale=1, color=['green', 'black'])
    plt.quiver(*origin, Q.T[:,0], Q.T[:,1], angles='xy', scale_units='xy', scale=1, color=['r', 'b'])
    plt.axis('square')
    plt.xlim([-2, 4])
    plt.ylim([-2, 4])
    plt.grid(True)
    plt.xlabel(r"$x$")
    plt.ylabel(r"$y$")
    plt.show()
```

### Decomposición $QR$ reducida


```python
def QR(A, modified=True):
    m, n = A.shape
    Q = np.zeros((m, n))
    R = np.zeros((n, n))
    # Gram-Schmidt orthogonalization
    for j in range(n):
        y = A[:, j]
        if modified: 
            for i in range(j):
                R[i, j] = np.dot(Q[:, i], y)
                y = y - R[i, j] * Q[:, i]
        else:
            for i in range(j):
                R[i, j] = np.dot(Q[:, i], A[:, j])
                y = y - R[i, j] * Q[:, i]
        R[j, j] = np.linalg.norm(y)
        Q[:, j] = y / R[j, j]
    return Q, R
```

### Ejemplo en $\mathbb{R}^2$

Para considerar el manejo de arreglos de ```numpy```, recordar que el procedimiento lo realizamos por vectores columnas, es decir,
\begin{equation}
    A=
    \left[
        \begin{array}{c|c}
            & \\ 
            \mathbf{a}_1 & \mathbf{a}_2 \\
            & \\
        \end{array}
    \right].
\end{equation}


```python
A = np.array([[2., 1.], [1., 3.]])
```


```python
A
```




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



#### Versión implementada


```python
Q, R = QR(A)
```


```python
np.dot(Q, R)
```




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



#### Versión ```numpy```


```python
Q2, R2 = np.linalg.qr(A)
```


```python
np.dot(Q2, R2)
```




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




```python
print("Implementado")
print("Q:\n", Q)
print("R:\n", R)
print()
print("Numpy")
print("Q:\n", Q2)
print("R:\n", R2)
```

    Implementado
    Q:
     [[ 0.89442719 -0.4472136 ]
     [ 0.4472136   0.89442719]]
    R:
     [[2.23606798 2.23606798]
     [0.         2.23606798]]
    
    Numpy
    Q:
     [[-0.89442719 -0.4472136 ]
     [-0.4472136   0.89442719]]
    R:
     [[-2.23606798 -2.23606798]
     [ 0.          2.23606798]]


La implementación de ```numpy``` no realiza la consideración de $r_{j,j} > 0$.

### ¿Son ortogonales?


```python
np.dot(Q[:,0], Q[:, 1]), np.dot(Q2[:,0], Q2[:, 1])
```




    (0.0, -5.551115123125783e-17)



Al parecer sí :)

### ¿Se cumple $Q^{-1}=Q^*$?


```python
np.linalg.inv(Q), Q.T
```




    (array([[ 0.89442719,  0.4472136 ],
            [-0.4472136 ,  0.89442719]]),
     array([[ 0.89442719,  0.4472136 ],
            [-0.4472136 ,  0.89442719]]))




```python
np.linalg.inv(Q2), Q2.T
```




    (array([[-0.89442719, -0.4472136 ],
            [-0.4472136 ,  0.89442719]]),
     array([[-0.89442719, -0.4472136 ],
            [-0.4472136 ,  0.89442719]]))



Al parecer también.

### Visualmente


```python
plotVector(A, Q)
```

Los vectores verde y negro corresponde a los vectores columnas de $A$. Los vectores rojo y azul corresponde a los vectores ortonormales encontrados por $QR$.

# Ajuste de curvas

## Resolución de mínimos cuadrados utilizando $QR$


```python
def QRLS(A, b):
    Q, R = QR(A) # QR reduced
    x = spla.solve_triangular(R, np.dot(Q.T, b)) # Solve Rx=Q^*b
    return x
```


```python
def plot(t, y, tt, yy, yyy):
    plt.figure(figsize=(12, 6))
    plt.plot(t, y, 'r.')
    plt.plot(tt, yy, label="Normal Equations")
    plt.plot(tt, yyy, label="QR")
    plt.grid(True)
    plt.legend()
    plt.xlabel(r'$t$')
    plt.ylabel(r'$y$')
    plt.show()
```

## Ecuaciones normales


```python
lm = lambda c_1, c_2, t: c_1 + c_2 * t
qm = lambda c_1, c_2, c_3, t: c_1 + c_2 * t + c_3 * t ** 2
em = lambda c_1, c_2, t: c_1 * np.exp(c_2 * t)
sm = lambda c_1, c_2, c_3, c_4, t: c_1 + c_2 * np.cos(2 * np.pi * t) + c_3 * np.sin(2 * np.pi * t) + c_4 * np.cos(4 * np.pi * t)
pm = lambda c_1, c_2, t: c_1 * t ** c_2
```


```python
# A^*A x = A^* b
def normalEquations(A, b):
    return np.linalg.solve(np.dot(A.T, A), np.dot(A.T, b))
```

#### Modelo lineal


```python
lm = lambda c_1, c_2, t: c_1 + c_2 * t
```


```python
n = 200
t_a, t_b = 0, 3
t = np.linspace(t_a, t_b, n)
tt = np.linspace(t_a, t_b, 3 * n)
```


```python
y_l = lm(1, 2, t) + np.random.normal(1.5, .5, n) # Data
A_l = np.ones((n, 2))
A_l[:,1] = t
b_l = y_l
```


```python
p_l = normalEquations(A_l, b_l)
```


```python
p_l_qr = QRLS(A_l, b_l)
```


```python
p_l, p_l_qr
```




    (array([2.61872557, 1.92680618]), array([2.61872557, 1.92680618]))




```python
plot(t, y_l, tt, lm(*p_l, tt), lm(*p_l_qr, tt))
```

Podemos ver para ambos métodos obtenemos los mismos coeficientes del modelo lineal.

## Cuadrático


```python
y_q = qm(1, 2, 3, t) + np.random.normal(5, 3, n)
A_q = np.ones((n, 3))
A_q[:, 1] = t
A_q[:, 2] = t ** 2
b_q = y_q
```


```python
p_l = normalEquations(A_q, b_q)
```


```python
p_l_qr = QRLS(A_q, b_q)
```


```python
plot(t, y_q, tt, qm(*p_l, tt), qm(*p_l_qr, tt))
```

## Periódico


```python
y_s = sm(1, 2, 2, 1, t) + np.random.normal(10, 2, n)
A_s = np.ones((n, 4))
A_s[:, 1] = np.cos(2 * np.pi * t)
A_s[:, 2] = np.sin(2 * np.pi * t)
A_s[:, 3] = np.cos(4 * np.pi * t)
b_s = y_s
```


```python
p_l = normalEquations(A_s, b_s)
```


```python
p_l_qr = QRLS(A_s, b_s)
```


```python
plot(t, y_s, tt, sm(*p_l, tt), sm(*p_l_qr, tt))
```

## Exponencial


```python
y_e = em(1, 2, t) * np.random.normal(10, 2, n)
A_e = np.ones((n, 2))
A_e[:, 1] = t
b_e = np.log(y_e)
```


```python
p_l = normalEquations(A_e, b_e)
```


```python
p_l_qr = QRLS(A_e, b_e)
```


```python
plot(t, y_e, tt,  em(np.exp(p_l[0]), p_l[1], tt),  em(np.exp(p_l_qr[0]), p_l_qr[1], tt))
```

## Comentarios

Podemos obtener los mismos resultados al utilizar **Ecuaciones Normales** que al usar la **Descomposición $QR$**. En teoría hay un método que requiere menos computación, ¿puede mostrarlo experimentalmente?

## Ejemplo 4.16 del libro guía p. 219
### ¿Por qué Gram-Schmidt modificado sobre el clásico?


```python
d = 1e-10
Ab = np.array([
    [1, 1, 1],
    [d, 0, 0],
    [0, d, 0],
    [0, 0, d]
])
```

### Clásico


```python
Qbc, Rbc = QR(Ab, False)
```


```python
print("Q:\n", Qbc)
print("R:\n", Rbc)
```

    Q:
     [[ 1.00000000e+00  0.00000000e+00  0.00000000e+00]
     [ 1.00000000e-10 -7.07106781e-01 -7.07106781e-01]
     [ 0.00000000e+00  7.07106781e-01  0.00000000e+00]
     [ 0.00000000e+00  0.00000000e+00  7.07106781e-01]]
    R:
     [[1.00000000e+00 1.00000000e+00 1.00000000e+00]
     [0.00000000e+00 1.41421356e-10 0.00000000e+00]
     [0.00000000e+00 0.00000000e+00 1.41421356e-10]]


¿$Q^*Q = I$?


```python
QTQ_c = np.dot(Qbc.T, Qbc)
QTQ_c
```




    array([[ 1.00000000e+00, -7.07106781e-11, -7.07106781e-11],
           [-7.07106781e-11,  1.00000000e+00,  5.00000000e-01],
           [-7.07106781e-11,  5.00000000e-01,  1.00000000e+00]])



¿$\|Q^*Q - I_n\|_F$?


```python
print(np.linalg.norm(QTQ_c - np.eye(3)))
```

    0.7071067811865477


### Modificado


```python
Qbm, Rbm = QR(Ab)
```


```python
print("Q:\n", Qbm)
print("R:\n", Rbm)
```

    Q:
     [[ 1.00000000e+00  0.00000000e+00  0.00000000e+00]
     [ 1.00000000e-10 -7.07106781e-01 -4.08248290e-01]
     [ 0.00000000e+00  7.07106781e-01 -4.08248290e-01]
     [ 0.00000000e+00  0.00000000e+00  8.16496581e-01]]
    R:
     [[1.00000000e+00 1.00000000e+00 1.00000000e+00]
     [0.00000000e+00 1.41421356e-10 7.07106781e-11]
     [0.00000000e+00 0.00000000e+00 1.22474487e-10]]


¿$Q^*Q = I$?


```python
QTQ_m = np.dot(Qbm.T, Qbm)
QTQ_m
```




    array([[ 1.00000000e+00, -7.07106781e-11, -4.08248290e-11],
           [-7.07106781e-11,  1.00000000e+00, -1.66533454e-16],
           [-4.08248290e-11, -1.66533454e-16,  1.00000000e+00]])



¿$\|Q^*Q - I_n\|_F$?


```python
print(np.linalg.norm(QTQ_m - np.eye(3)))
```

    1.1547005383859233e-10


Experimentalmente podemos ver que en este ejemplo los vectores de $Q$ obtenidos con el algoritmo clásico no son ortonormales. Notamos que con el algoritmo modificado no obtenemos $I_n$ exactamente, pero el resultado es mejor.

b0ed908258313ea62dbc4766de20356077bf155c
ipynb
Jupyter Notebook
material/06_minimos_cuadrados/qr.ipynb
fitoxgs2/INF-285
f7d547c2480aea7eb1242b873ab1135a0873834f
BSD-3-Clause
2020-04-24T01:25:24.000Z
2021-05-10T01:08:37.000Z
material/06_minimos_cuadrados/qr.ipynb
fitoxgs2/INF-285
f7d547c2480aea7eb1242b873ab1135a0873834f
BSD-3-Clause
2020-04-22T00:57:29.000Z
2020-04-23T23:59:28.000Z
material/06_minimos_cuadrados/qr.ipynb
fitoxgs2/INF-285
f7d547c2480aea7eb1242b873ab1135a0873834f
BSD-3-Clause
2020-04-20T01:15:09.000Z
2020-07-27T22:53:59.000Z
Qwen/Qwen-72B
1. YES
2. YES