4
$\begingroup$

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

Improve this question
$\endgroup$
5
  • $\begingroup$ how large might be $N$? is $\mathbf{A}$ (even though its elements depend only on $i,j$) sparse? $\endgroup$
    – Anton Menshov ♦
    Commented Nov 20, 2018 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
    Commented Nov 20, 2018 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
    Commented Nov 21, 2018 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
    Commented Nov 21, 2018 at 14:50
  • $\begingroup$ Looks like a job for the fast multipole method? $\endgroup$
    – Federico Poloni
    Commented Nov 24, 2018 at 23:34

1 Answer 1

Reset to default
3
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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