Is there any library or routine for high-performance matrix-matrix product, where the matrix elements are computed on-the-fly using a given function of $i$ and $j$?

More specifically, in the problem I am currently facing, I have to compute a matrix-matrix product

$$ \mathbf{y} = \mathbf{A} \mathbf{b} $$

where $\mathbf{A}$ is a symmetric matrix $N \times N$, $\mathbf{y}$ and $\mathbf{b}$ are $N\times 3$.

Using LAPACK or BLAS routines requires the formation of the matrix $\mathbf{A}$. The matrix elements, however, are a function that depends only on the matrix indices $i,j$.

  • $\begingroup$ how large might be $N$? is $\mathbf{A}$ (even though its elements depend only on $i,j$) sparse? $\endgroup$ – Anton Menshov Nov 20 '18 at 17:28
  • 2
    $\begingroup$ How big is your $N$? Is $A$ sparse or dense? How large a submatrix of $A$ can you work with at one time in main memory? How simple is the function that computes the entry of $A$ for a given $i$ and $j$? $\endgroup$ – Brian Borchers Nov 20 '18 at 20:01
  • $\begingroup$ @BrianBorchers The matrix is dense. N will be at the order of 1,000 to 5,000. I do not form the matrix $\mathbf{A}$; all elements are computed on-the-fly; thus I am not restricted to the size of the submatrix that can fit in main memory. The function is $\mathbf{A}(i,j) = 1/(1+\|x_i - x_j\|^2)$ where $x$ is a known vector (size $N\times 1$). This is not the only function, but an example of a function I am currently using. $\endgroup$ – fcdimitr Nov 21 '18 at 12:29
  • $\begingroup$ You haven't given any context to explain why it is necessary to avoid forming the matrix $A$? Your matrix is quite small and could easily fit into memory. Since $A_{i,j}$ is quite expensive to compute it will be fastest to just form $A$ and use optimized BLAS routines to do the multiplication. $\endgroup$ – Brian Borchers Nov 21 '18 at 14:50
  • $\begingroup$ Looks like a job for the fast multipole method? $\endgroup$ – Federico Poloni Nov 24 '18 at 23:34

If the matrix is small enough to fit into memory, then there is of course no cost associated with actually forming the elements: You will have to compute the elements at least once anyway to perform the matrix-vector product, so you might as well store the values you computed in memory and then run an optimized matrix-vector product function over it.

But there are of course cases where you genuinely can't or don't want to compute and store the matrix entries. An example is if $A=BC$ where $B$ and $C$ are other matrices and where it is cheaper to compute $y=Ab$ via $B(Cb)$ than $(BC)b$. Other cases are where $A$ does not fit into memory.

In these cases, you need a library that allows you to implement a matrix via its action, i.e., its matrix-vector product. One example is PETSc where you can create a MatShell object -- a matrix for which you need to provide a pointer to a function that can perform matrix-vector products. Most other linear algebra libraries have similar functionality -- I'm pretty sure Trilinos can do this, and deal.II (disclaimer: that's a project I'm associated with) does this all over the place by creating classes that provide a matrix-vector multiplication member function.


Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.