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

Question

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

Answer 1

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.

Answer 2

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}\\ 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$.

Answer 3

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.

Answer 4

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