0
$\begingroup$

I am solving the following matrix equation for $\mathbf{x}$:

$$(J^{\mathbf{T}}J)\mathbf{x}=J^{\mathbf{T}}\mathbf{r}$$

  • $J$ is $m\times n$ matrix
  • $\mathbf{x}$ is vector of size $n$
  • $\mathbf{r}$ is vector of size $m$

Furthermore, $m\gg n$.

It is sometimes the case that $n$ is quite large (e.g. 1000) and $J$ contains values like 500. Standard multiplication is then impossible to do even with double-precision numbers as the values exceed $10^{128}$.

I found that scaling $J$ by a nonzero value has predictable effect on the result (it has the same effect as multiplying both sides of the above equation by scalars), so I am thinking about computing a scale factor so that the above equation can be solved without overflows. However, too small scale factor can also endanger accuracy of the matrix multiplication.

Furthermore, $J$ can contain both very small and very large values.

So the question is how to choose some scaling factor $s$ so that we can substitute $J$ by $sJ$ and avoid too big numbers and also keep high accuracy of the multiplication.

My ideas:

  • make average/median of elements in $J$ equal to e.g. $10^3$
  • scale $J$ according to sum of its biggest column
  • find some "preconditioner" matrix $A$ and substitute $J$ for $AJ$ ... the preconditioned solution can be obtained using inverse transform, but $A$ can get very large as it has to be $m\times m$

Any ideas?

$\endgroup$
  • $\begingroup$ Double precision numbers have a range of $10^{\pm 308}$, far exceeding the value you wrote ($10^{128}$). Could you clarify the maximum your values can actually reach? $\endgroup$ – Doug Lipinski Mar 6 '15 at 17:14
  • $\begingroup$ Okay so its larger than that... there are 1000 numbers in range 600-1000 multiplied and added together. $\endgroup$ – Libor Mar 6 '15 at 17:41
  • $\begingroup$ The question can be simplified to: How to solve the equation given that $m$ is large and $J$ contains large values. $\endgroup$ – Libor Mar 6 '15 at 17:43
  • $\begingroup$ @Doug Okay the problem was that $J$ contained NaN (not-a-number) values. I have not noticed this. This explains the result - it was not casued by an overflow. $\endgroup$ – Libor Mar 6 '15 at 18:24
  • $\begingroup$ If you compute the dot product of two vectors of length 1000 containing numbers each as large as 1000, you get 1000*1000*1000 = 10^9. Where does $10^{128}$ comes from? $\endgroup$ – Federico Poloni Mar 6 '15 at 19:48
4
$\begingroup$

This isn't a direct answer to your question but rather an alternative approach.

It looks like you are solving a least squares (LS) problem using the normal equations. The normal equations are known to be a poor way to solve least squares problems because $J^T J$ is often very ill-conditioned.

So why don't you simply solve the LS problem directly ($J x = r$) using, for example, a QR decomposition? The LAPACK routine, DGELS, will do that for you.

http://www.netlib.org/lapack/lug/node27.html

There are interfaces to this routine in Octave, MATLAB, Python, and other scripting languages so you don't need to use Fortran/C/C++ if you'd prefer not to.

If you are interested in reading more about using QR factorizations to solve LS problems, this chapter in Cleve Moler's book provides a very readable introduction:

http://www.mathworks.com/moler/leastsquares.pdf

$\endgroup$
  • $\begingroup$ Yes. But I am not using simple normal equations, but rather SVD of $J^{\mathbf{T}}J$ inside Levenberg-Marquadt step (this approximates Hessian, which is very convenient and the matrices are small for large datasets, it also converges much faster that 1st order methods). But I wanted to keep the question simple and focus on the core problem. I cannot use LAPACK as I am working on a commercial project that cannot have any extra dependencies. It is also routine that works perfectly unless big numbers appear. $\endgroup$ – Libor Mar 6 '15 at 16:48
  • 1
    $\begingroup$ Same thing, you can easily derive the SVD of $J^TJ$ from that of $J$. There should be no need to compute the product explicitly. $\endgroup$ – Federico Poloni Mar 6 '15 at 19:45
  • $\begingroup$ In fact you don't want to compute $J^{T}J$ because that squares the condition number. Just use the AVD of $J$ directly (or more cheaply use a QR factorization) $\endgroup$ – Brian Borchers Mar 6 '15 at 20:57

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.