Re: matrix library ?

From: lums_at_[hidden]
Date: 2001-03-17 08:49:02

--- In boost_at_y..., "David Abrahams" <abrahams_at_m...> wrote:
> 
> ----- Original Message -----
> From: <lums_at_l...>
> To: <boost_at_y...>
> Sent: Friday, March 16, 2001 9:00 AM
> Subject: [boost] Re: matrix library ?
> 
> 
> > --- In boost_at_y..., "David Abrahams" <abrahams_at_m...> wrote:
> > > A might not be an explicitly represented
> > > > matrix of any kind -- it could simply be a function that 
returns
> > the
> > > > result A*x without using a stored matrix (e.g., a matrix-free
> > > > operator).
> > >
> > > Not much code in apply() in that case ;-)
> >
> > Haha!  Actually there would be something.  The typical way of 
doing a
> > matrix-free linear operator application is to use our old friend
> > Taylor's expansion
> >
> >   A*y ~= [ f(x + d*y) - f(x) ] / dy
> >
> > So that A is the Jacobian of the nonlinear function f.  To 
represent
> > A requires storing f, x, and the result of f(x) and then doing the
> > above computation.  There will also be some computing to size dy 
just
> > right so as to make it as small as possible but not too small,
> > trading off mathematical accuracy for floating point accuracy.
> 
> Forgive my denseness, but the math chops are rusty!
> a. That's not a Taylor series expansion the way I'm used to seeing 
it. Is
> there some trivial derivation I'm missing?

Taylor is in there -- f(x+dy) ~= f(x) + f'(x)*dy -- in this case, f 
is a function from R^n to R^n so f' is a matrix (the Jacobian of f).  
We can rewrite

f(x+dy) ~= f(x) + A*dy

where A is the Jacobian -- rearranging gives what I had before 
(modulo a typo):
  
  A*y ~= [ f(x+dy) - f(x) ] / d

for some small (but not too small) d.

What this is doing is the following.  If you have a nonlinear 
function and you want to compute the matrix-vector product based on 
the Jacobian of this function times a vector, you do not need to 
construct any kind of matrix to do this.  You just need to evaluate 
the non-linear function itself at some perturbed location (x+dy).

> b. Which is the linear operator?

A, the Jacobian of the non-linear function f.

 
> and,
> 
> 3. You'd want to use more straightforward code if you had actual 
matrices,
> no?

Of course.  This was just to illustrate that there are cases (and 
this technique is actually used) where one might have a linear 
operator that is not a stored matrix of values.

Quick trivia moment.  The etymology for 'matrix' is 'womb.'  This 
makes much sense in the context of the movie with that title.  I am 
not sure how it applies to linear algebra.....