# Minimizing Cost Functions using Iterative Least Squares

I am currently trying to use iterative least squares to solve a system, $y = Hx + v$ where $y$ is a vector of observations, $H$ is the design matrix, and $v$ is the observation error.

From my understanding, iterative least squares tries to minimize a cost $J$ where $J = \min (Hx-y)^2$.

My pseudocode for this is

while ||x^k-x^{k-1}||
H = h(x^{k-1})
r = f(x^{k-1}) - y
dx = pinv(H)*r
x^{k} = x^{k-1} + dx
end


where $H$ is the Jacobian and the output of $f(x)$ is the expected observations.

Can $r$ be called the cost of the system because it is the residual/difference between $f(x)$ at estimate $x^{k-1}$ and $y$? Should it be squared in an implementation of iterative least squares?

• Since $r$ is apparently a vector, you probably meant to ask if $||r||$, the norm of the residual, can be called "the cost of the system". As with $J$, cost is usually a scalar. – hardmath Aug 10 '15 at 11:51

The linear least squares problem $\min \|Hx-y\|^2$ is equivalent to solving $H^THx = H^Ty$, known as the normal equation. $H^T$ denotes the transpose, or adjoint matrix. There are many ways to solve this, depending on your problem, but the conjugate gradient method is a good general method.
Note that the pseudo-inverse (pinv) function you're using will calculate $(H^TH)^{-1}H^T$, so
This means your iterative method would have converged in the first iteration. If performance is not an issue and $H$ is reasonably small, that is the easiest way to solve your problem, but for larger $H$ iterative methods will often be significantly faster. This is particularly the case if $H$ is a sparse matrix, or the matrix-vector product $Hr$ can be calculated via a fast method such as the Fourier transform.
As hardmath mentioned, $r$ is the residual vector and $J=\|r\|^2$.
• As a follow-up, would iterating $(H'WH)^{-1}H'Wy$ be effectively similar since an inverse is being calculated? Would the conjugate gradient method be a good approach for sparse $H$ matrices? – d.mc2 Aug 10 '15 at 13:13
• Yes, you would just be solving the weighted least squares, instead of least squares. Conjugate gradient is an excellent general method when $H$ is sparse, particularly if a good preconditioner can be found. – LKlevin Aug 11 '15 at 12:05