# Avoid arithmetic overflow in matrix multiplication

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?

• 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? Mar 6 '15 at 17:14
• Okay so its larger than that... there are 1000 numbers in range 600-1000 multiplied and added together. Mar 6 '15 at 17:41
• The question can be simplified to: How to solve the equation given that $m$ is large and $J$ contains large values. Mar 6 '15 at 17:43
• @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. Mar 6 '15 at 18:24
• 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? Mar 6 '15 at 19:48

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

• 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. Mar 6 '15 at 16:48
• 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. Mar 6 '15 at 19:45
• 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) Mar 6 '15 at 20:57