# Fast multiplication of constant symmetric positive-definite matrix and vector.

Consider the matrix $H=H^T$, $H>0$, $H \in R^{n \times n}$, and the vector $v \in R^n$. In a numerical algorithm, I need to compute the product $b = Hv$. Right now I am following the naive approach: $b_i = \sum_{j=1}^{n} h_{ij} v_j, i=1,...,n$. Is there a faster way to compute this product? $H$ is non-sparse and constant (i.e. eigenvectors, eigenvalues, etc. of $H$ are available).

-
Thank you everyone. In some of the computing platforms on which the algorithm will run (e.g. microcontroller without FPU), there is no BLAS implementation –  user19781 Dec 7 '11 at 16:00

If I understand the question right, by "constant" it is meant that $H$ is a fixed, but arbitrary positive definite matrix.

In general, I don't think that you can compute the matrix-vector product $Hv$ faster than $O(n^2)$. But if $H$ has structure (Toeplitz, Circulant, Strictly diagonally dominant, etc.), or if you are willing to settle for approximate answers, you can compute the product faster ("randomized linear algebra" is a good search term).

A very naive approach is the following: Suppose that $v$ is some vector, and instead of using $H$, you use $H_k$, the top-$k$ rank approximation to $H$. Then, $H_kv$ can be computed in time $O(nk)$. The error of this computation is $\|Hv-H_kv\| \le \|H-H_k\|\times\|v\|$, which is fairly easy to characterize.

-
$H$ is fixed and has no special structure. Thanks for the search term. –  user19781 Dec 7 '11 at 15:56

In practice no, unless $H$ has some additional structure that you haven't told us about (low rank blocks, displacement structure, or similar stuff).

There are some results on faster-than-$O(n^3)$ matrix multiplication (by the way, the exponent has recently been lowered to $O(n^{2.373})$), but as far as I know none of them is actually viable in practice.

Grab the fastest BLAS implementation that you can find for your machine, and use it.

-
Federico: isn't the OP asking for BLAS2? –  Suvrit Dec 7 '11 at 11:35
@Suvrit: I have always heard and used $\text{BLAS}=\{\text{BLAS level 1}, \text{BLAS level 2}, \text{BLAS level 3}\}$. Or am I misunderstanding what you're saying? –  Federico Poloni Dec 7 '11 at 11:40
Actually, as far as I know the Strassen algorithm IS viable in practice for reasonably large sized matrices (how large is "reasonably large" depends on how vectorized your FPU is, historically 100x100 was big enough). This will help the OP, since he can perform $n$ of his multiplications in time $(O(n^{2.7}),$ so $O(n^{1.7})$ per multiplication. –  Igor Rivin Dec 7 '11 at 11:46
@Federico: I meant that the OP is asking for a matrix-vector multiplication, not matrix-matrix. But as Igor also notes, as do you, if the OP requires $Hv$ for several $v$'s then, yes, it makes sense to start thinking about matrix-mult. algos. –  Suvrit Dec 7 '11 at 12:48

In fact, $O(n^2)$ arithmetic operations is not avoidable, in general. A general result by Winograd (see "Algebraic complexity theory" by Buergisser, Clausen, and Shokrollahi, Sect. 13.2) shows that the for a generic $m\times n$ matrix (all the entries are algebraically independent) and a generic vector one will need $O(mn)$ multiplications and additions (there is a precise formula too, but $O(mn)$ would do for our purposes). PSD matrices are symmetric, so directly it won't apply t $H$. However, you can cut $H$ into blocks: $H=\begin{pmatrix} H_{11}&H_{12}&bsol;\ H_{12}^\top&H_{22}\end{pmatrix}$. Then $H_{12}$ is generic, in general, and, writing $v$ as a block vector $v=(v^1,v^2)$ one has $b=(b^1,b^2)=(H_{11}v^1+H_{12}v^2, H_{12}^\top v^1+H_{22}v^2)$. Now Winograd's result is applicable to $H_{12}$, so you still get $O(n^2)$, as you cannot avoid computing $H_{12}v^2$ when you try to get $b$.

-

In practice, on typical desktop computers and server class machines using the x86-64 architecture, matrix-vector multiplication is limited more by memory bandwidth than floating point operations. This happens because there are no opportunities in matrix-vector multiplication for bringing data in from memory to cache and reusing it before flushing it out of the cache. Getting a matrix-vector multiply to run at full memory bandwidth is generally quite easy.

As others have pointed out, you might be able to go faster if the matrix has specialized structure (e.g. if it's a Toeplitz matrix you can do the multiplications in O(n*log(n)) time.)

-
You can use the symmetry of your matrix to cut the memory traffic in half, but it takes careful coding to get that efficiency in practice. –  Brian Borchers Dec 7 '11 at 15:55