How do we compute a matrix-vector product without using the standard matrix multiplication?

Question

I am a bit confused right now. I am taking a class on numerical linear algebra and as I have understood sometimes one doesn't compute the Matrix-Vector product $Av$ normally but uses other techniques to get the result.

What are these techniques and are they so much faster than the "standard way"?

Additional Info:

I think that this description of the problem was a bit misleading. What I actually wanted to know was how one computes $Av$ without actually forming A(I realized my confusion by reading the answers of @Amit Hochman and @Laurent90).

Answer 1

Many linear algebra problems, such as solving a linear equation, finding eigenvalues, etc., can be solved iteratively provided you can compute the matrix-vector product. For example, in the power iteration, if you compute $v_{n+1} = Av_n$, starting from a random vector, you will generally arrive at the eigenvector corresponding to the largest eigenvalue. In many practical situations, $A$ will be very large, perhaps too large to fit in memory, but nonetheless, its product by a vector can be computed efficiently, without ever forming the matrix.

Say, for example, $A = D + uu ^T$, where $D$ is diagonal and $u$ is a column vector, of length $n$. Then you compute the matrix-vector product as $Av = Dv + u(u^Tv)$, which can be done in O(n) space and time complexity, without ever forming the matrix. Such sparse + low-rank matrices occur often in practice. Another example is that instead of $uu^T$ you have a discrete-Fourier transform matrix, and you can compute its matrix-vector product using the FFT. There are numerous algorithms along these lines, the most well-known of which is the fast multipole method (https://en.m.wikipedia.org/wiki/Fast_multipole_method).

Answer 2

It is possible that your $A$ is the Jacobian of a residual vector function $f$ which you need to drive to zero to solve your problem. It can also be the Jacobian corresponding to space or space-time discretisation of a PDE.

For instance for the 1D heat transfer $\partial_t T = D \partial_{xx} T$ on a uniform grid, applying a finite difference scheme in space could be, for the $i$-th space point:

$$\partial_t T_i = \dfrac{D}{\Delta x^2} (T_{i+1} -2 T_{i} + T_{i-1})$$

In matrix form, you get: $$\partial_t \mathbf{T} = \dfrac{D}{\Delta x^2}L \mathbf{T} = A \mathbf{T}$$

With $L$ the classical Laplacian matrix.

Usually, for more complex problems, you do not form $A$, sometimes you may not even know how it's formed. Instead you write a function which directly computes the time derivatives at each point as per the first discrete equation. Computing $A\mathbf{T}$ in that case is equivalent to calling this residual function with T as input, which may be both simpler and faster.

That's for instance how JNFK approaches work.

I can come back later to give a simple Python approach of both approaches if that helps clarifying.

EDIT: I was typically thinking of the 1D heat transfer. Here is a small Python example that shows to ways of semi-discretising the PDE with finite differences, with Neumann conditions on the sides:

import numpy as np
from scipy.sparse import diags

N = 50 # number of mesh points
L=1 # domain length
dx = 1/(N-1) # space step
T0=np.linspace(0,1,N)**2 # initial solution

# 1st approach, form the Laplacian matrix
A = diags([1,-2,1],[1,0,-1], shape=(N,N)).toarray() # not good, drops the sparse character
A[0,0]=-1 # Neumann condition
A[-1,-1]=-1
A = A/dx**2
def dTdt_matrix(T):
  return A.dot(T)

# 2nd approach, directly compute the terms pointvby point
def dTdt(T):
  dTdt = np.zeros((T.size,))
  dTdt[1:-1] = (1/dx**2) * ( T[:-2] - 2*T[1:-1] + T[2:] )
  dTdt[0]  = (1/dx**2) * ( T[1]  - T[0]  )
  dTdt[-1] = (1/dx**2) * ( T[-2] - T[-1] )
  return dTdt

assert np.allclose( dTdt(T0), dTdt_matrix(T0) )

The second approach is faster, and is very easy to code. It is also more error proof in my experience (generalising to 2D for instance). This is even more true if the time derivative is a nonlinear function. Note however that there are instances where it may be needed anyway, or at least approximate the Jacobian of $f$ (here $A$), for instance for implicit time integration. Jacobian-free Newton-Krylov approaches do not compute $A$, but only matrix vector products which can be obtained with the model function directly. $A$ or some approximation of $A$ may however be needed for preconditioning.

Answer 3

The most common place I can think of where you don't actually compute a matrix-vector product is in solving the linear problem $$A x = b$$ Naively, $$ x = A^{-1} b $$ which requires you to first compute the inverse of $A$ before performing a matrix-vector product. The matrix-vector product isn't necessarily expensive (though it's not free, either), however computing $A^{-1}$ is very expensive compared to other methods which can directly solve for $x$, for example Gaussian elimination or LU decomposition coupled with a triangular matrix solver.