# 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?

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)$$.

• Is $n=m$ in your case? Otherwise $L^k$ is not defined for $k\ge 2$. – iljusch Dec 2 at 0:39
• Good point. Yes, let's stick with $n=m$. (I will edit the question soon). – user124784 Dec 2 at 4:24
• Quick follow-up, we could alternatively think about finding $(LL^T)^k$ for $k>=1$ in the even case (defined similarly for the odd case). But the square case is interesting enough for me. – user124784 Dec 2 at 4:32
• Here is evidence that this is not easy. Let $W =\left[ \begin{smallmatrix} 1&1 \\ 0 &1 \end{smallmatrix} \right]$. Then $U L^2 V^T = \frac{1}{\sqrt{5}} \left[ \begin{smallmatrix} 3&4 \\ 1&3 \end{smallmatrix} \right]$. So we have to leave the world of integer operations. On the other hand, a square root is not that bad ... – DES-SupportsMonicaAndTransfolk Dec 2 at 15:33

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 over-arching problem I see to this is whether or not this actually runs down your computational complexity. This might be an interesting read for you.

• This is a nice observation! I will check out the paper you provided (before the bounty expires) and see if this will solve my problem. At the very least, it has given me a new angle to look at this. Thanks! – user124784 Dec 3 at 7:49
• The matrix square root observation is interesting, but unfortunately I don't think this really solves the problem at its core. Specifically, this method requires the square root to be instantiated explicitly which seems to be at least an O(n^3), requiring matrix-matrix multiplication. I haven't been able to find a way to compute A^{1/2}x more cheaply than computing A^{1/2}. Ultimately, we want to take advantage of the fact that we only need matrix-vector products. – user124784 Dec 4 at 17:07
• @user124784 That's fair enough. It's quite difficult to devise good, computationally efficient matrix expressions for algorithms. Best of luck with your research! – Mnifldz Dec 5 at 1:21

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!

• Interesting approach! This method would certainly be more efficient than the SVD, giving O(kdmn) in the end. If good approximations can be found for relatively small d (which seems to be the case), then this could work well for me. I'll try this one out soon! – user124784 Dec 4 at 17:17
• You're welcome. Also, $O(dmn)$, not $O(kdmn)$. The exponent $k$ comes in only in how it effects the choice of $d$. – DES-SupportsMonicaAndTransfolk Dec 4 at 20:59
• Ah yes, you're right! I've awarded the bounty to you as, at the very least, this has given me some promizing things to push on. I'll select as an answer once I've implemented and confirmed that it does actually do what it says on the tin. Thanks! – user124784 Dec 5 at 1:34