4
$\begingroup$

We have a large sparse symmetric positive-definite matrix $A \in \mathbb R^{N \times N}$ and a vector $x \in \mathbb R^N$. How do I practically compute the inner product $x^T A x$ when the matrix $A$ is ill-conditioned?

In my applications, this is a stiffness or mass matrix of a finite element code.

There is an obvious algorithm that first computes $y = A x$ and then the dot product $x^T y$. However, this product occasionally happens to be negative in my computations. That cannot be in exact arithmetic but it happens in the presence of rounding errors.

I am wondering how that can be mitigated.

Improve this question
$\endgroup$
5
  • 3
    $\begingroup$ Have you tried using Kahan or Neumaier summations in place of regular summations when computing the matrix-vector products, and then again performing the dot product? $\endgroup$
    – helloworld922
    Commented Mar 4 at 21:32
  • 1
    $\begingroup$ Is it small enough to perform SVD? If you compute coordinates $y$ of $x$ in the principal basis of $A$, then the task becomes to compute $\|y\|^2$. If $A$ is ill-conditioned, some entries of $y$ will be huge and others tiny, so you have to add them in correct order to avoid the "small+large" computation. $\endgroup$
    – Yaroslav Bulatov
    Commented Mar 4 at 21:35
  • 1
    $\begingroup$ How do you know for sure that your matrix is positive definite? Where does it come from? Do you have an explicit expression as $A = BB^*$ available maybe? $\endgroup$
    – Federico Poloni
    Commented Mar 4 at 22:51
  • 2
    $\begingroup$ If you can decompose $A$ into $R^*R$ then $x^*R^*Rx = \|Rx\|^2$ which will never be negative. $\endgroup$
    – lightxbulb
    Commented Mar 5 at 12:25
  • $\begingroup$ I suspect your $A$ or $x$ has a wide range of values, a stable way might be to normalize $x$ or $A$ and then reorder the indices according to the magnitude of your input values, which can reduce the rounding a little bit. Otherwise, when you take the summation, try the reduction algorithm, in practice, it is not much worse than the Kahan summation if $x$ is only at million level. $\endgroup$
    – Yimin
    Commented Mar 5 at 19:42

1 Answer 1

Reset to default
8
$\begingroup$

You don't need a temporary vector. Instead, you loop over the elements of the matrix $(i,j)$ and update a single counter:

  sum = 0
  for ( (i,j) in matrix entries )
    sum += A(i,j) * x(i) * x(j)

This algorithm works for both dense and sparse matrices.

In practice, this appears to be stable to within round-off for the matrices that appear in the finite element method. In 25 years of working with and on finite element software, I have not encountered the need to use things such as Kahan summation for the sum above. Instead, if I've encountered issues where this sum should have been but wasn't positive, it was invariably a bug in assembling the matrix.

Improve this answer
$\endgroup$
6
  • 1
    $\begingroup$ I suspect this might be a real issue in cases in which the mass matrix is only positive semidefinite, however. $\endgroup$
    – Federico Poloni
    Commented Mar 5 at 18:00
  • $\begingroup$ when $A(i,j)$ are of comparable sizes, and $x$ does not have a huge range of values, the summation should be fine since the rounding is almost a random variable and sometimes there is no need to round off. I guess the OP encounters a case that the rounding is a huge issue. $\endgroup$
    – Yimin
    Commented Mar 5 at 19:35
  • $\begingroup$ It is not that uncommon for this to happen, e.g. with very ill-conditioned matrices. Incidentally I was just reading Gould, Hribar, and Nocedal's "ON THE SOLUTION OF EQUALITY CONSTRAINED QUADRATIC PROGRAMMING PROBLEMS ARISING IN OPTIMIZATION" paper, and there they explicitly mention such a case: "In both cases the CG iteration was terminated when $r^Tg$ became negative, which indicates that severe errors have occurred since $r^Tg$ must be positive." $\endgroup$
    – lightxbulb
    Commented Mar 5 at 21:11
  • $\begingroup$ @FedericoPoloni The Laplace matrix can become semidefinite, but the mass matrix better does not. That would imply that you have cells of zero area. $\endgroup$
    – Wolfgang Bangerth
    Commented Mar 5 at 22:08
  • $\begingroup$ @Yimin If $A$ is the mass matrix, then $x$ is a solution vector. Solutions can vary by a few orders of magnitude, but I cannot imagine a case where the solution of a PDE would vary by such a factor that you get into trouble with summation. $\endgroup$
    – Wolfgang Bangerth
    Commented Mar 5 at 22:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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