# Fast matrix multiplication with matrix elements computed on-the-fly (without forming the matrix)

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

• how large might be $N$? is $\mathbf{A}$ (even though its elements depend only on $i,j$) sparse? – Anton Menshov Nov 20 '18 at 17:28
• 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$? – Brian Borchers Nov 20 '18 at 20:01
• @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. – fcdimitr Nov 21 '18 at 12:29
• 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. – Brian Borchers Nov 21 '18 at 14:50
• Looks like a job for the fast multipole method? – Federico Poloni Nov 24 '18 at 23:34

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.