4
$\begingroup$

Currently, I encountered a problem of approximating the following series:

$$ (I-X)^{-\frac{1}{2}}=I+\frac{1}{2}X+\frac{1\cdot3}{2\cdot4}X^{2}+\frac{1\cdot3\cdot5}{2\cdot4\cdot6}X^{3}+\ldots $$

where $X$ is a diagonalizable matrix and the largest (in absolute value) eigenvalue is less than $1-\frac{1}{\kappa},\kappa>1$. What I am looking for is something like the following

$$ (1-X)^{-1}=I+X+X^{2}+X^{3}+\dots\approx\prod_{k=0}^{d}(1+X^{2^{k}}) $$

With $O(\log(\kappa\log(1/\epsilon)))$ multiplications and summations, we can achieve accuracy of $O(\epsilon)$.

We can assume the following representations are cheap to get (we assume that their computational cost is $1$):

  1. $1+aX^{n},a\in(-1,1),n\in[-1,0,1,\dots)$

  2. Multiplication between any two cheap representations.

  3. Summation between any two cheap representations.

So how can we approximate $(I-X)^{-\frac{1}{2}}$ efficiently?

My initial idea is

$$ (I-X)^{-\frac{1}{2}}=I+\frac{1}{2}X(1+\frac{3}{4}X(1+\dots(1+\frac{2n-1}{2n}X(1-X)^{-1}))\dots) $$

But it is not efficient enough.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

Not exactly the same thing, but you can use Newton's method easily in your framework. As you can see for instance here, Newton's method for the inverse square root can be computed with additions and multiplications only, and the same algebraic expression that works for scalars can be used for diagonalizable matrices without modifications (proof: in the eigenvector basis, doing the iteration on the matrix level is equivalent to running $n$ scalar iterations inside a diagonal matrix).

It converges quadratically, so it should have the convergence rate that you need.

By the way, your initial idea is known in numerical analysis as Horner's scheme for polynomial evaluation.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thank you very much for your answer. It partially solved my problem. Although the quadratic convergence of Newton's method is exactly what I am looking for, but expanding the Newton's update ($d$ steps,e.g., from $y_0=I$ to $y_1$ to $y_d$ ) will result in $2^d$ times matrix operations. The reason I don't want to store $y_1$ to $y_d$ is that the matrix multiplication is usually as expensive as calculating its SVD and sorry for not making this clear in the question. Your answer definitely gives me something to work on. $\endgroup$
    – cdh
    Commented Apr 4, 2014 at 19:58
  • $\begingroup$ @DehuaCheng I see. Don't you have the same problem in your example with $(I-X)^{-1}$, then? If you expand out all the products, it becomes $2^d$ operations again. $\endgroup$
    – Federico Poloni
    Commented Apr 5, 2014 at 6:57
  • $\begingroup$ Yes, you are right. In that inverse case, by applying sparsifier on $X^2,X^4,\dots,X^{2^d}$, we can actually afford these multiplications, so it takes $d$ multiplications to get $X^2,X^4,\dots,X^{2^d}$. In the inverse square root case, I am now trying to do the similar thing by sparsifying $S_n$ (assuming the update rule is $y_{n+1} = y_n S_n$). There are still difficulties, but it looks promising. $\endgroup$
    – cdh
    Commented Apr 6, 2014 at 7:34

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .