[PYTHON] Least squares method (triangular matrix calculation)

Hello. The least squares method is $ J (\ boldsymbol {x}) \ triangleq \ left (\ boldsymbol {y} --A \ boldsymbol {x} \ right) ^ \ textrm {T} \ left (\ boldsymbol {y} --A \ boldsymbol {x} \ right) Calculates $ \ boldsymbol {x} _ {\ min} $ with $ as the least. It is a normal way to go through a triangular matrix for good calculation efficiency. The details are the following relational expressions, but to summarize the calculation,

Prepare the matrix $ \ tilde {A} \ triangleq \ begin {bmatrix} A & \ boldsymbol {y} \ end {bmatrix} $ (Ay in the example below) and then (square) up Triangular matrix $ \ tilde {R} \ triangleq \ begin {bmatrix} R & \ tilde {\ boldsymbol {y}} \\ \ boldsymbol {0} ^ \ textrm {T} & \ alpha \ end {bmatrix} $ ( In the example below, find `` `Ry```) and
Use this to solve $ R , \ boldsymbol {x} = \ tilde {\ boldsymbol {y}} $ to find $ \ boldsymbol {x} _ {\ min} $.

An example Python implementation of this is [^ 1],

def lsqsolv(Ay):
    n = Ay.shape[1] - 1
    Ry = triangulateupper(Ay)
    return numpy.linalg.solve(Ry[0:n,0:n], Ry[0:n,n]) # x_min

def triangulateupper(A):
    return numpy.linalg.cholesky(numpy.dot(A.transpose(),A)).transpose()

[^ 1]: The triangular matrix `triangulateupper ()` in this example uses the Cholesky decomposition, which emphasizes computational efficiency, but in many cases, the QR decomposition, which emphasizes numerical stability, is recommended.

To summarize the relational expressions ($ \ tilde {Q} $, $ Q $ are orthogonal matrices),

\begin{align*}
J (\boldsymbol{x}) &\triangleq \left(\boldsymbol{y} - A \boldsymbol{x} \right)^\textrm{T} \left(\boldsymbol{y} - A \boldsymbol{x} \right) \\
&= \begin{bmatrix} \boldsymbol{x}\\ -1\end{bmatrix}^\textrm{T} \tilde{A}^\textrm{T} \, \tilde{A} \,\begin{bmatrix}  \boldsymbol{x}\\ -1\end{bmatrix} \\
&= \begin{bmatrix} \boldsymbol{x}\\ -1\end{bmatrix}^\textrm{T} \tilde{R}^\textrm{T} \, \tilde{R} \,\begin{bmatrix}  \boldsymbol{x}\\ -1\end{bmatrix} \\
A &= Q R \\
\tilde{A} &= \tilde{Q} \tilde{R} \\
\tilde{R}^\textrm{T} \, \tilde{R} &= \tilde{A}^\textrm{T} \tilde{A} \\
\tilde{A} &\triangleq \begin{bmatrix} A & \boldsymbol{y} \end{bmatrix} \\
\tilde{R} &\triangleq \begin{bmatrix} R  & \tilde{\boldsymbol{y}} \\ \boldsymbol{0}^\textrm{T} & \alpha \end{bmatrix} \\
\tilde{\boldsymbol{y}} &\triangleq Q^\textrm{T} \boldsymbol{y} \\
\alpha &\triangleq \left( \min_{\boldsymbol{x}} J (\boldsymbol{x}) \right)^{1/2} \\
 \\
\end{align*}