# Solve $\left( \sum_i J_i J_i^T \right) x = \sum_i J_i \cdot \alpha_i$ efficiently

Is it possible to solve for $x$

$$\left( \sum_i^N J_i J_i^T \right) x = \sum_i^N J_i \alpha_i$$ where $\alpha_i \in \mathbb{R}$ is a scalar and $J_i \in \mathbb{R}^6$ is a vector

in such a way that I do not have to do a summation over all $\mathbb{R}^{6 \times 6}$ matrices ($J_i J_i^T \in \mathbb{R}^{6 \times 6}$).

### What is the reason?

• $N \approx 10^5$ is quite large.
• On a GPU this involves 2 sum reductions.
• The left-hand side involves a $\mathbb{R}^{6 \times 6}$ sum reduction - the right-hand side a $\mathbb{R}^{6}$ sum reduction: $$\left( \sum_i^N \mathbb{R}^{6 \times 6}_i \right) x = \sum_i^N \mathbb{R}_i^{6} \mathbb{R}_i$$
• Reducing $\mathbb{R}^{6 \times 6}$ matrices $N$ times drains the GPU bandwidth.

### What is my goal?

• Doing a sum reduction with a lower dimensionality. For example $$x = \sum_i^N \mathbb{R}^{6}$$

### Important notes

• The GPU is memory bound but not compute bound.
• Consider $\left( \sum_i^N J_i J_i^T \right)$ as positive semi-definite.

Thanks

• As in your previous post, it is not clear that $\left(\sum_iJ_iJ_i^T\right)^{-1}$ even exists. Depending on the number and values of the $J_i$ there could be infinitely many solutions, exactly one solution, or no solution at all.
– cthl
Commented Jan 12, 2017 at 1:39

If you can find a vector $X$ that is perpendicular to all but one of the $J_i$'s, then you can left multiply by $X^T$ to get $$\sum_{i} X^T J_i J_i^T x = \sum_{i} X^T J_i g,$$ which will simplify to $$X^T J_k J_k^T x = X^T J_k g$$ where we assume that $X^T J_k \neq 0$. This simplifies to $$J_k^T x = g$$ which, as in your previous question, has infinitely many solutions.

• I updated the question. Please revise. Thanks :) Commented May 31, 2017 at 14:46

$\sum_i J_i J_i^T$ is a positive semidefinite matrix, positive definite if the inverse exists. The Cholesky decomposition is then then both efficient and numerically stable.

• Thanks, but the problem is not really solving the system of equation. It's more the reduction step. Commented May 31, 2017 at 14:45

I don't understand your rationale exactly, but what about conjugate gradient? This will solve your problem exactly in exactly $6$ steps (and $6$ gradient evaluations/ matrix $A=\sum J_iJ_i^T$ vector multiplies). You don't need to actually compute $A$ itself, you only need to compute $Ax$, for vectors, which can be done by computing the individual $J_i\times(J_i\cdot x)$ terms in the sum.

• That sounds promising, can you elaborate a little please? Commented Jun 6, 2017 at 10:05
• en.wikipedia.org/wiki/Conjugate_gradient_method All you need to do is evaluate gradients and other much smaller calculations to do conjugate gradient. There is a theorem somewhere that says that conjugate gradient will converge to the exact solution in exactly $N$ steps, where $N$ is the number of dimensions. IN your case $N=6$, so this will be rather fast. Commented Jun 6, 2017 at 16:06
• What kind of elaboration do you need? Commented Jun 6, 2017 at 16:38