
Boost : 
From: lums_at_[hidden]
Date: 20010317 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 matrixfree
> > > > operator).
> > >
> > > Not much code in apply() in that case ;)
> >
> > Haha! Actually there would be something. The typical way of
doing a
> > matrixfree 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 matrixvector 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 nonlinear function itself at some perturbed location (x+dy).
> b. Which is the linear operator?
A, the Jacobian of the nonlinear 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.....
Boost list run by bdawes at acm.org, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk