2
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – hardmath Aug 10 '15 at 11:51
1
$\begingroup$

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

x=pinv(H)*y

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$.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – d.mc2 Aug 10 '15 at 13:13
  • $\begingroup$ 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. $\endgroup$ – LKlevin Aug 11 '15 at 12:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.