# Can these matrices be multiplied in $\mathcal O(n^2)$ time?

Consider a real-valued orthogonal matrix $$Q$$ and a sequence of diagonal matrices $$\{D_m\}_{m=1}^\infty$$. All entries of $$Q$$ are real and the entries of each $$D_n$$ are real and positive. What is the cost of the following multiplication(s)?

$$\Sigma_m = Q D_m Q^\top$$

Is there any way this can be accomplished in $$\mathcal O(n^2)$$ time? I'm interested in the amortized cost, meaning that I'm okay with a $$\mathcal O(n^3)$$ pre-processing step if it leads to (eventual) $$\mathcal O(n^2)$$ multiplications.

• This doesn't answer the question, but there are fast algorithms to multiply a matrix with its transpose. So you could factor $D$ (ending up with complex numbers) and then compute $Q \sqrt{D}$ in $n^2$, and then the full product in time $2/(2^\omega-3) n^\omega$. Jul 8, 2020 at 3:32
• @tch that's a very helpful observation - thank you. Jul 8, 2020 at 18:04
• What is a (stochastic) diagonal matrix? I understand that to be a matrix whose rows (or columns) sum to 1, which would imply an identity matrix. Jul 8, 2020 at 18:08
• @AlexR. I meant that the matrix is itself is a stochastic process. I forgot that row/column stochastic has a very specific meaning in matrix algebra. I have edited this out of the question as I don't think its relevant to the question. All that matters is that we don't know $D_{n+1}$ until after we have computed $\Sigma_n$. Jul 8, 2020 at 18:11
• Out of curiosity, what's your goal with getting the original matrix? I'm just wondering if maybe you can do what you truly want without reconstructing the original matrix? Jul 8, 2020 at 22:02

I'm not aware of any algorithm for this, although perhaps someone knows of one. Nothing I write below uses the fact that $$Q$$ is orthogonal, so perhaps that can be exploited in some way.

In terms of theoretical asymptotic complexity, $$AA^T$$ can be computed slightly faster than an arbitrary matrix product can be. Specifically, faster by a factor $$2/(2^\omega-3)$$ where $$\omega$$ is currently about 2.37.

However, since your matrices are of size $$n=1000$$, the constants suppressed by the big-O, and more importantly, the implementation of the matrix product algorithms will come into play. Effectively using a library such as numpy which takes advantage e of CPU caching, factorization, etc. will likely have a very large impact on the runtime.

## symmetry

If we're assuming we are going to pay $$n^3$$, then we should focus on the constant in front. An obvious optimization is that $$QD_mQ^T$$ is symmetric, so only the upper (or lower) triangular part needs to be computed. This can be done directly bu noting that, $$[\Sigma_m]_{i,j} = Q_{[:,i]}^T D_m Q_{[:,j]} = \sum_{k=1}^{n} [D_m]_{k,k} Q_{i,k} Q_{k,j}$$ where $$Q_{[:,j]}$$ is the $$j$$-th column of $$Q$$. Naively computing $$\Sigma_m$$ this way would require $$3n-1$$ flops, and need to be done for $$n(n+1)/2$$ entries. If $$m=1,\ldots, M$$, then this is a total cost of $$M(3n-1)n(n+1)/2\approx (3/2)Mn^3$$.

## preprocessing

However, notice that the only dependence of this sum on $$D_m$$ is the term $$[D_{m}]_{k,k}$$. Thus, the product $$Q_{i,k}Q_{k,j}$$ can be precomputed. If we define a new vector $$q^{(i,j)} = [Q_{i,1}Q_{1,j}, Q_{i,2}Q_{2,j}, \ldots, Q_{i,n}Q_{n,j} ]^T$$ then we have $$[\Sigma_{m}]_{i,j} = (q^{(i,j)})^T d_m$$, where $$d_m$$ is the vector with entries equal to those of the diagonal of $$D_m$$.

By symmetry, $$q^{(i,j)} = q^{(j,i)}$$, so we can preprocess our data set by computing $$q^{(i,j)}$$, $$i=1,2,\ldots, n$$, $$j=1,2,\ldots, i$$, and storing the $$n(n+1)/2$$ many vectors. The cost of computing each $$q^{(i,j)}$$ is $$n$$, so the total cost is $$n^2(n+1)/2 \approx n^3/2$$ (also this much storage is required).

## vectorization

Now the cost of computing $$[\Sigma_{m}]_{i,j}$$ is the cost of a dot product: $$2n-1$$ which must be done $$n(n+1)/2$$ times for each $$\Sigma_m$$. This gives a total cost of $$n^2(n+1)/2 + M(2n-1)n(n+1)/2 \approx n^3/2 + M n^3$$, so if $$M$$ is very large we improve by a factor of roughly 1.5.

But, perhaps even more important is that all of these products can be vectorized. Specifically, we can store the $$q^{(i,j)}$$s in a $$n(n+1)/2\times n$$ matrix, and all the $$d_m$$s in a $$n\times M$$ matrix. The last step would be to take the data out of this matrix product and put it into a useable form.

What is faster in practice will probably depend a lot on what libraries you're able to use, and how they manage memory, etc.

• Great answer - thanks. Unfortunately, I can't use vectorization because the matrix $D_{m+1}$ isn't given until after we have "evaluated" the product at time $m$. Jul 10, 2020 at 4:51