 # Boost :

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