0
$\begingroup$

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

$\endgroup$
7
  • $\begingroup$ Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking. $\endgroup$
    – Community Bot
    Jul 28, 2022 at 16:28
  • $\begingroup$ Google Strassen or Winograd algorithm $\endgroup$
    – davidhigh
    Jul 28, 2022 at 16:45
  • $\begingroup$ Also, you might want to explain what you mean by "doesn't compute the matrix-vector product $Av$ normally". I can think of 20 different ways of computing these products, depending on context and what information I have, and none of them strike me as "normal" :-) $\endgroup$ Jul 28, 2022 at 22:47
  • $\begingroup$ Any investigation of algorithms must begin with a specification of the inputs and the output. Supposing for the sake of discussion that matrix $A$ and vector $v$ of compatible sizes are inputs and the output sought is product vector $Av$, there's nothing inherently wrong with looping through the rows of $A$ to get the column $Av$ one entry at a time. Improvement on that "standard definition" approach will necessarily require some information about the inputs $A$ and $v$. $\endgroup$
    – hardmath
    Jul 29, 2022 at 2:54
  • 1
    $\begingroup$ @davidhigh None of the algorithms you cite applies to matrix-vector products; just to matrix-matrix ones (and I don't think they provide any advantage when the last dimension is 1). $\endgroup$ Jul 29, 2022 at 11:24

3 Answers 3

5
$\begingroup$

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

$\endgroup$
4
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ I think this is exactly what I was looking for. You are writing: "...you do not form A, sometimes you may not even know how it's formed". How can I imagine computing "Av" when I don't even know how to form A? Do you maybe mean that the Matrix A could be so big, that there is no way to save it as a MAtrix in the computer and therefore we need other techniques? $\endgroup$
    – Josh.K
    Jul 29, 2022 at 12:52
  • 1
    $\begingroup$ @Josh.K I guess you can always form $A$ in theory, but it may be costly or tedious in some instances. I've updated my anwser with the heat equation as an example. $\endgroup$
    – Laurent90
    Jul 29, 2022 at 18:55
2
$\begingroup$

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.

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