Given $W=ULV^T$ and a vector $\mathbf{x}$, can we compute $UL^kV^T\mathbf{x}$ without doing the SVD, for any integer k?

Question

Consider a matrix $W \in \mathbb{R}^{n\times m}$ with corresponding singular value decomposition, $W = ULV^T$, and a vector $\mathbf{x} \in \mathbb{R}^{m}$. Is it possible to compute matrix-vector-products of the form $U L^k V^T \mathbf{x}$ without explicitly computing the SVD?

If $k$ is odd, this is easy. For example, $W W^{T} W\mathbf{x} = UL^3V^T\mathbf{x}$. But I can't see an obvious way to solve this for even powers.

For context, I want to avoid the SVD as computing the matrix vector products will typically be much cheaper (e.g. for $k \ll n)$. Ideally, a working algorithm would have complexity $O(kmn)$ instead of $O(mn^2 + nm^2)$.

Answer 1

I highly doubt there is an exact formula. But, if you want a good approximation, here is an idea: Find a polynomial $f(x) = \sum f_j x^j$ of degree $d$ which is a good approximation of $x^{(k-1)/2}$. (More on this later.) Then $\sum f_j (WW^T)^j W \vec{x} = U \left( \sum f_j L^{2j+1} \right) V^T \vec{x}$. For each singular value $\lambda$, $\sum f_j \lambda^{2j+1}$ will be $\lambda f(\lambda^2) \approx \lambda^k$. If we organize the evaluation as in Horner's method as $W( f_0 \vec{x} + \cdots + W^T W(f_{d-2} \vec{x} + W^T W (f_{d-1} \vec{x} + W^T W f_d \vec{x})) \cdots )$, it will only take $O(dmn)$ steps.

Now, how to approximate $x^{k-1/2}$ by polynomials? Minor observations first: (1) We can precompute our $f(x)$ for many values of $k$, so this needn't count against our computation time if we expect $k$ to stay bounded. (2) A variant of this might be to precompute a good approximation $g(x)$ for $x^{(k-1)/2}$ on $[0,1]$, then compute a bound $M$ for the squares of the singular values, such as $\mathrm{Tr}(W^T W)$, and use $f(x) = M^{(k-1)/2} g(x/M)$. (3) We could just take an approximation $h(x)$ to $\sqrt{x}$ and multiply by $x^{k/2-1}$, or raise it to the $2k-1$ power, but neither of those is likely to be the best option; it would probably be better to optimize for each $k$ separately.

With that out of the way, how well can polynomials approximate $\sqrt{x}$? A numerical analyst would have better strategies, but what I did was to expand $\sqrt{x}$ in the Laguerre polynomials. I chose Laguerre polynomials because they are natural for functions defined on $[0, \infty)$, but I am not an expert and, if this is important, I would suggest trying several families of orthogonal polynomials to see which is better.

The cubic approximation is $$\frac{\sqrt{\pi }}{192} \left(x^3-15 x^2+90 x+30\right)$$ and does an okay job on the range $[0,4]$ and not so well on $[4,9]$ (Laguerre in blue, $\sqrt{x}$ in red).

The degree $10$ approximation is $\sqrt{\pi}$ times

(167610643200 + 1676106432000 x - 1257079824000 x^2 + 670442572800 x^3 - 209513304000 x^4 + 39109150080 x^5 - 4444221600 x^6 + 306979200 x^7 - 12471030 x^8 + 271700 x^9 - 2431 x^10)/1902536294400

and looks pretty good!

Answer 2

I think an important observation to make here, especially since the SVD is derived from the spectral theorem, is that the matrices $WW^T$ and $W^TW$ are symmetric positive semi-definite matrices. As mentioned by iljusch in the comments, we will stick to the case that $n = m$. In particular we find that

$$ WW^T \;\; =\;\; UL^2U^T \hspace{3pc} W^TW \;\; =\;\; VL^2V^T. $$

Because they are positive semi-definite we can in fact find their square-roots since all of the singular values will be nonnegative:

$$ \sqrt{WW^T} \;\; =\;\; ULU^T \hspace{3pc} \sqrt{W^TW} \;\; =\;\; VLV^T. $$

Therefore in computing even powers we find that square-roots can serve as a key to performing your computations without actually performing SVD: \begin{eqnarray*} \sqrt{WW^T}W & = & UL^2V^T \\ \sqrt{WW^T}(WW^TW) & = & \left (WW^T\right )^{3/2}W \\ & = & UL^4V^T. \\ \end{eqnarray*}

Looks as if your general formula for positive integers is given by: $$ UL^kV^T \;\; =\;\; \left (WW^T\right )^{\frac{k-1}{2}}W. $$

This can just as well be formulated as $$ UL^kV^T \;\; =\;\; W \left (W^TW\right )^{\frac{k-1}{2}}. $$

The one remaining problem I see is whether or not this ultimately lowers your computational complexity. This might be an interesting read for you.