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
> > 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
> > 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
> > 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
> > 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
> 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.
> 3. You'd want to use more straightforward code if you had actual
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.....
Boost list run by bdawes at acm.org, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk