mm-184
===
Subject: Re: Kalman filter help> Can you suggest a
book that is not too heavy on the math? I am mostly> looking
for implementation examples. I have looked at enough>
derivations and would rather look at real world examples.>>Try
R.M. Rogers, >Applied Mathematics in Integrated Navigation
Systems>, >>published by AIAA. It is all about using Kalman
'lters to combine >>navigation solutions from multiple
sensors.>>Lars>>
===
>>Subject: Re: Kalman 'lter help>>Try this
URL>>
http://ourworld.compuserve.com/homepages/PDJoseph/kalman.htm>>
I bought the lessons from Peter Joseph and the seemed to be
quite clear.> Can you suggest a book that is not too heavy
on the math? I am mostly>> looking for implementation
examples. I have looked at enough>> derivations and would
rather look at real world examples.>> Try R.M. Rogers,
>Applied Mathematics in Integrated Navigation Systems>, >
published by AIAA. It is all about using Kalman 'lters to
combine > navigation solutions from multiple sensors.>>
Lars>>
===
>>Subject: How to choose 3 unity vectors from 1000
unity vectors, the volume of the parallelepiped formed by the
dimensional unity vectors, one can form a>>parallelepiped
using arbitrary 3 of them.>>I want to choose the 3 unity
vectors from the 1000 unity vectors, the>>volume of the
parallelepiped spanned by these 3 unity vectors should>>be the
largest amongst all the parallelepiped formed by any possible
3>>unity vectors in the 1000 unity vectors.>>I cannot choose
the speci'ed 3 unity vectors by calculate the volumes>>of all
parallelepiped formed by all possible combination of 3
unity>>vectors from the 1000 unity vectors, because the
combination number is>>too large: about 1.6e8.>>Is there any
for your attention.>>Kenki>>The Hong Kong Polytechnic
University.>>
===
>>Subject: Re: How to choose 3 unity vectors
from 1000 unity vectors, the volume of the parallelepiped
All:>>Given 1000 three dimensional unity vectors, one can
form a>parallelepiped using arbitrary 3 of them.>>I want
to choose the 3 unity vectors from the 1000 unity vectors,
the>volume of the parallelepiped spanned by these 3 unity
vectors should>be the largest amongst all the parallelepiped
formed by any possible 3>unity vectors in the 1000 unity
vectors.>>I cannot choose the speci'ed 3 unity vectors by
calculate the volumes>of all parallelepiped formed by all
possible combination of 3 unity>vectors from the 1000 unity
vectors, because the combination number is>too large: about
1.6e8.>>Why do you believe that 1.6e8 is too large ?>>What
computer are you using ? Mine solves the problem in 16 seconds
cpu time.>>
===
>>Subject: Re: How to choose 3 unity vectors from
1000 unity vectors, the volume of the parallelepiped formed by
the 3 unity vectors are the largest?>> kenkicn@yahoo.com
unity vectors, one can form a>>parallelepiped using arbitrary
3 of them.>I want to choose the 3 unity vectors from the
1000 unity vectors, the>>volume of the parallelepiped spanned
by these 3 unity vectors should>>be the largest amongst all
the parallelepiped formed by any possible 3>>unity vectors in
the 1000 unity vectors.>I cannot choose the speci'ed 3 unity
vectors by calculate the volumes>>of all parallelepiped formed
by all possible combination of 3 unity>>vectors from the 1000
unity vectors, because the combination number is>>too large:
about 1.6e8.>Is there any easier way to choose the speci'ed
Hong Kong Polytechnic University.>>for any vector in the set:
take this as the 'rst. perform a householder
>>orthogonalization step transfering it to the 'rst unit
vector and apply this >>transformation to the other vectors
too. >>then, select the set of vectors of maximal>>length with
respect to >>component 2 and 3 under the transformed ones. for
those in this set, repeat by >>transforming >>the lower
subvector to a multiple (less than one in abs value) of the
>>second unit vector i.e. the form (*,*,0)>>and apply the
transformation also to the remaining 998, select one of
largest >>third component and take those three. proof:
maximizes the possible value of the>>absolute value of the
determinant = volume, since transformation applied does >>not
change the absolute value of the determinant of the selected 3
by 3 system.>>by observing the fact that the diagonal elements
in this qr-decomposition >>will decrease in absolute value you
need not do all the work and will be able>>to skip a lot, but
with respect to the 'rst vector i see no chance to avoid
>>testing them all (you could have the unit matrix within your
system somewhere>>but the remaining vectors being nearly
linearly dependent from one of the unit>>vector for
example)>>hth>>peter>>
===
>>Subject: Re: SVD algoritm
wanted...>>X-Beer: Yes please>>X-Woobinda: Oh yes>>X-Kebab:
9Z-b=R4zSD6](xNB_#I[:xDc3)_[S4>G]i[&O:GAo'WLOf?@v;lA%=u1TXx:{
g6$]p|KO]kMPMl1;EQ41g+dzM>nb|YP}o@8U(,.n4{fDyxrw#n7ek2U-!$b^u*
,e55g4F=!v?GJ#%:su7E&E>> Laz <usenet@b-tarts.co.uk>
either written in>>C or which can easily be called from C?>>I
have been using the svd algorithm from Numerical Recipes in C
but now,>>as I need to distribute some code, I need to replace
it with a free>>alternative.>>Is the best way to use one of
the subroutines from LAPACK? If this is>>the case, can anyone
point me in the direction of some simple examples>>where
people have called these subroutines from C? I have never
called>>Fortran from within C so I'm 'nding my 
feet here a
bit!>>Are there any other alternatives which I could
consider?>> clapack in http://www.netlib.org> or>
c/meschach in the same source have svd in c.>>The latter looks
just what I'm after: a pretty much self-contained
function.>>This will eventually be for distribution so I 
don't
===
>>Subject:
Re: SVD algoritm wanted...> Laz <usenet@b-tarts.co.uk>
algorithm either written in>>C or which can easily be called
from C?>>I have been using the svd algorithm from
Numerical Recipes in C but now,>>as I need to distribute
some code, I need to replace it with a free>>alternative.>>
>Is the best way to use one of the subroutines from
LAPACK? If this is>>the case, can anyone point me in the
direction of some simple examples>>where people have called
these subroutines from C? I have never called>>Fortran from
within C so I'm 'nding my feet here a bit!>>Are 
there
any other alternatives which I could consider?>> clapack in
http://www.netlib.org>> or>> c/meschach in the same source
have svd in c.>>The latter looks just what I'm after: a
pretty much self-contained function.>This will eventually be
for distribution so I don't want lots of
svd() routine has a bug that manifests itself>>in the sorting
of the results ( vectors and values ). It>>can be fatal -- it
crashed ( some of ) my programs -- depending>>on the platform,
etc.>>Basically, the Meschach non-recursive sorting routine
taken>>from Sedgewick c1990 _Algorithms_ ( or one of the
>Algorithms>texts ) propagates an error. There is a problem
where the>>current singular value being sorted is compared to
a value>>held by one of the sentinels. A conditional statement
is>>missing that allows the singular value comparison
during>>sorting to run amok. It's an easy 'x. 
I'll append
some>>notes from my e-mail to Stewart from 1/97, a note to
myself>>and some data to illustrate the problem.>>AFAIK,
Sedgewick has corrected the omission in later texts.>>The
Meschach svd() routine has proven robust
otherwise.>>HTH.>>/david>>-------------------- %cut here%
------------------->>Here's where the problem occurs in
'xsvd():>> ...>> /* i = partition(d->ve,l,r) */>> ...>> { /*
inequalities are >backwards> for *>> while ( d->ve[++i] > v
)>> ;>> while ( d->ve[--j] < v)>> ;>> ...>>The data:>> 1
1732.8370 2156.4650>> 1 2142.9590 1034.9940>> 1 2132.1760
1056.9140>> 1 2201.4720 1061.8640>> 1 2119.4480 809.7800>> 1
1504.2650 1830.4890>> 1 1381.2280 2085.4010>> 1 1330.4930
2144.2670>>-------------------- %cut here%
------------------->>-------------------- %cut here%
------------------->>the internal non-recursive quicksort
needs>>an additional check (conditional) per>>the
sentinels.>>original.>> while ( d->ve[++i] > v ) ;>> while (
d->ve[--j] < v ) ;>>one solution. we negate the values to be
consistent>>with sedgewick's examples.>> while ( -d->ve[++i] 
<
-v ) if ( i == r ) break ;>> while ( -d->ve[--j] > -v ) if ( j
== l ) break ;>>this is slight overkill. instead,>> while (
-d->ve[++i] < -v ) ; /* d --> v ; OK */>> while ( -d->ve[--j]
> -v ) if ( j == l ) break ;>>that is, d->ve[++i] increasing
approaches the partitioning>>element v since v =
d->ve[d->dim-1] . in the case>>of d->ve[--j] the loop may
continue below the de'ned>>array d to values of j that are
negative. boo.>>-------------------- %cut here%
------------------->>/david>>
===
>>Subject: Re: mrqmin
parameter constraints>> jjjj7777@hotmail.com (Jinesh)
all you people. So>>I will eliminate the simpli'cation and
explain the entire problem.>>Here's the rig:>I have an array
of data points (x) from my experiment. I would like to>>'t a
bi-exponential equation to the data points.>Equation: y =
A*exp(t/TA) + B*exp(t/TB)>Independent params are A, TA, B,
TB>>t is an array for the times for which I have the data
points.>>x is array of experimental data points>>A & B are in
range [0 2000]>>TA & TB are in range [0 500]>Objective:
Obtain values for A, TA, B, TB such that (y - x) is
the>>minimum. In other words, the best 't.>Like I said,
Matlab constrained non-linear minimization
function,>>lsqcurve't, gives me the correct results. mrqmin()
gives me incorrect>>results since atleast one the parameters
is negative, thereby>>rendering the entire solution
meaningless.>Any help, gurus?>>if you already have a
solution, why ask?>>if you had looked at the page given
already two times then you would have found >>also the special
exponential 't code exp't (in matlab). if you 
need >>a solution
different from matlab and cannot use f77, then gauss't 
>>would
be the way to go (as already said). otherwise, since your
>>problem has even more structure, namely separability (you
can express >>A and B as functions of TA and TB using the
necessary optimality >>conditions) there is even a better way,
namely using dqed from Hanson >>and krogh, which allows also
for bounds and uses this
structure.>>hth>>peter>
===
>>Subject: Re: mrqmin parameter
constraints> if you already have a solution, why ask?> if
you had looked at the page given already two times then you
would have found > also the special exponential 't code
exp't (in matlab). if you need > a solution different from
matlab and cannot use f77, then gauss't > would be the way
to go (as already said). otherwise, since your > problem has
even more structure, namely separability (you can express > A
and B as functions of TA and TB using the necessary optimality
> conditions) there is even a better way, namely using dqed
from Hanson > and krogh, which allows also for bounds and
uses this structure.> hth> peter>>Two reasons for asking.
The MATLAB lsqcurve't program does not>>calculate the
covariance matrix and I need to know the covariance of>>the
'tted parameters. Additionally, I want a C implementation of
the>>non-linear optimization function. Matlab was just for the
prototype>>testing. Matlab compiler does generate C binaries
but the compiler is>>centuries away from being optimized. I
cannot use that.>>I only learned today that mixing Fortan and
everyone for the suggestions.>>
===
>>Subject: sinc(x+y)
relations for sinc(x+y) like>>there are for sines and cosines
(i.e. sin(x+y) =>>sin(x)cos(y)+cos(x)sin(y) )? I don't think
so becasue of how the sinc>>function is de'ned (sin(x)/x) - 
or
it's not obvious if there are, but>>thought I would ask
===
>>Subject: Re:
angle-sum relations for sinc(x+y) like> there are for sines
and cosines (i.e. sin(x+y) => sin(x)cos(y)+cos(x)sin(y) )? I
don't think so becasue of how the sinc> function is 
de'ned
(sin(x)/x) - or it's not obvious if there are, but> thought
about , if>>Sinc(z,x) := xsin(x)/(z^2+x^2) and Cosc(z,x) :=
xcos(x)/(z^2+x^2),>>then>>sinc(x+y) =
(2/pi)Integrate[Sinc(z,x)Cosc(z,y)
+>>Sinc(z,y)Cosc(z,x),{z,0,In'nity}].>>-- >>HTH,>>Gerry
T.>>
===
suggestion. I played around with it a little, but 
wasn't>>able
to simplify things any further. But I did learn something new,
relations for sinc(x+y) like>> there are for sines and cosines
(i.e. sin(x+y) =>> sin(x)cos(y)+cos(x)sin(y) )? I don't 
think
so becasue of how the sinc>> function is de'ned (sin(x)/x) -
or it's not obvious if there are, but>> thought I would ask
if>> Sinc(z,x) := xsin(x)/(z^2+x^2) and Cosc(z,x) :=
xcos(x)/(z^2+x^2),>> then>> sinc(x+y) =
(2/pi)Integrate[Sinc(z,x)Cosc(z,y) +>
Sinc(z,y)Cosc(z,x),{z,0,In'nity}].>> -- > HTH,> Gerry
T.>>
===
>>Subject: Re: sinc(x+y) relations?>Pat>
for the suggestion. I played around with it a little, but
wasn't> able to simplify things any further. But I did learn
something new,>>since> I've never seen your relations
so thank you, it's the only stimulation 
I've>>had all day aside
from 2 coffees. Didn't Knuth remark that nowadays the>>only 
one
to get excited when a new identity is found is its discoverer.
See>>below for a slight simpli'cation that 
you've undoubtedly
already noticed .>>Gerry Thomas> <gfthomas@sympatico.ca>
angle-sum relations for sinc(x+y)>>like>> there are for sines
and cosines (i.e. sin(x+y) =>> sin(x)cos(y)+cos(x)sin(y) )? I
don't think so becasue of how the>>sinc>> function is 
de'ned
(sin(x)/x) - or it's not obvious if there are,>>but>> 
thought
about , if>> Sinc(z,x) := xsin(x)/(z^2+x^2) and Cosc(z,x) :=
xcos(x)/(z^2+x^2),>> then>> sinc(x+y) =
(2/pi)Integrate[Sinc(z,x)Cosc(z,y) +>>
Sinc(z,y)Cosc(z,x),{z,0,In'nity}].>>The foregoing can be
rewritten as>>sinc(x+y) = {sinc(x)cosc(y) +
sinc(y)cosc(x)}(2/pi)Integrate[(1 ->>z^2/(z^2+x^2))(1 -
z^2/(z^2+y^2)),{z,0,In'nity}].>>-- >>Ciao,>>Gerry
T.>>
===
>>Subject: Fast Fourier Transform>>Any one can help
me?>>Im a PhD student in Energy Engineerig, at this very
moment Im>>dealing with FFT for 'ltering data (water levels)
drawn from a power>>plant. For this Im using Matlab, Ive
obtained the power spectrum and>>the frequency of my data so
far, but I dont know how to get rid of>>the unwanted noise in
matlab. The amplitude goes from 0 to 0.025 and>>my frequency
varies from 0 to 250 Hz.>>What can I do to eliminate the
===
noise?.>>Subject: Re: Fast Fourier Transform>>
>Victor>
me?> Im a PhD student in Energy Engineerig, at this very
moment Im> dealing with FFT for 'ltering data (water levels)
drawn from a power> plant. For this Im using Matlab, Ive
obtained the power spectrum and> the frequency of my data so
far, but I dont know how to get rid of> the unwanted noise in
matlab. The amplitude goes from 0 to 0.025 and> my frequency
varies from 0 to 250 Hz.> What can I do to eliminate the
noise?.>>Apply some 'ltering method? If your noise is at high
frequencies and useful>>signal at lower frequencies simple
low-pass 'lter with appropriate>>bandwidth will do
it.>>
===
>>Subject: Modeling ground reaction forces (vehicle
I have yet to have an idea>>of what I'm actually up 
to.>>I'm in
the business of developing §ight simulators, and one part
of>>that is to model the forces and moments produced by the
landing gear>>when the aircraft lands and rolls on the ground.
>>The underlying physics is clearly the same for all
ground-bound vehicles>>like cars, so perhaps somebody from the
vehicle simulation community>>might be able to provide valuable
information.>>Generally, a landing gear is modeled by an
``open-loop'' approach:>>First, the forces and 
moments
produced by aerodynamics and propulsion>>are computed, then,
if some wheel touches the ground, the respective>>reaction
forces and moments. Finally, they are all summed up and
the>>state equation for velocities and positions of the
vehicle is>>integrated.>>Typically, gear models seem to
consist of one spring-damper strut>>obeying some quadratic
damping law, often one steerable or freely>>castoring wheel,
and tire and brake models of varying degrees of>>complexity
(I've even seen brake temperature modeled!).>>All these
approaches suffer from the notorious 
``drifting''
problem:>>When the vehicle is stopped on the ground (governed
by static friction),>>it slowly drifts and turns rather
randomly, while, in reality, it would>>remain stationary, at
most slightly moving in response to aerodynamic>>and
propulsion. More precisely, the points where the wheels touch
the>>ground move, while in reality they really remain
stationary as long as>>static friction governs (up to the
maximum static friction forces).>>The standard solution
employed for this problem seems to be 
to>>``freeze'' vehicle
position and orientation once relative speed (between>>ground
and vehicle) falls below a certain threshold.>>However, this
approach is quite ``un-physical'' and it is 
also dif'cult>>to
implement given a moving ground (imagine you model earth's
rotation,>>or an aircraft carrier rolling on the ocean waves).
And you also loose>>the above described swinging due to
aerodynamics or propulsion. This is>>quite real, a light
aircraft will typically just remain stationary but>>lower its
nose (pitch down) if full brakes and full power are
applied.>>In this situation, the vehicle reference point
_will_ (and should)>>slightly move.>>It would seem more
natural to me to freeze the points of ground contact.>>But,
you have more than one wheel and if you do this, you
actually>>restrict the way the vehicle can move:>> ___>> |>>
--------------------->> | |>> o o>>Considering this drawing,
if you 'x the points of ground contact of the>>two wheels
shown, you geometrically restrict the aircraft to move
only>>vertically (We assume a rigid vehicle apart from the
wheel struts). On>>the one hand, this is unrealistic, and on
the other hand this seems to>>me like asking for numerical
trouble (although I have to admit I have>>not yet written down
the equations). And this is only a
two-dimensional>>simplic'cation...>>This is assuming that the
struts expand and compress only vertically.>>Would it perhaps
help to allow also horizontal and lateral displacement>>with
some strong spring and damping constants?>>Well, these are, so
far, my thoughts on the subject. I'm sorry that>>this 
isn't
very concrete, but I want to make sure to head in the
right>>direction from the beginning when I tackle this
problem. I would very>>much appreciate discussion with experts
in this 'eld and any hints and>>references!>>Kind
trap, please see my homepage for my>>real address!>>--
>>Gerhard Wesp o o Tel.: +41 (0) 43 5347636>>Bachtobelstrasse
56 | http://www.cosy.sbg.ac.at/~gwesp/>>CH-8045 Zuerich
_/>>
===
>>Subject: Re: Modeling ground reaction forces (vehicle
very much appreciate discussion with experts in this 'eld
and> any hints and references!>>Well, I can start you off
with a reference: the last time I recall reading>>detailed
discussion of these sorts of problems, I was playing with
the>>Open Dynamics Engine (currently at
http://opende.sourceforge.net/ ). You>>may not want to
actually use their code, but you'll probably at least>>enjoy
the conceptual discussion in their user's guide.>>--->>Roy
Stogner>>
===
>>Subject: Re: Q. On Crout's Algorithm For Matrix
LU Decomposition>Fritz Foetzl> <fritz_foetzl@hotmail.com>
student, trying to write a software> library for handling
matrices that, among other things, 'nds> determinants
ef'ciently. The best way to do this is by LU> decomposition,
and the best way to do that (that I know of) is Crout's>
algorithm.>> I stumbled upon Crout's algorithm in _Numerical
Recipes In C: The Art> of Scienti'c Computing_, which has a
pretty good discussion. I> understand most of it.>> Where
I get stuck is the concept of the >pivoting> (Swapping
rows,> when appropriate. This results in a row-wise
permutation of the> the stability of Crout's method.>
Calculation of the subdiagonal> elements (alpha[i][j])
involves dividing by a pivot element. The> mathematical
discussion of the algorithm just uses the non-unity>
diagonal elements (beta[j][j]), without regard to pivoting. I
worked> through it by hand with a random 3x3 matrix, and it
seemed to work> just 'ne. However, the example C code goes
out of its way to choose> maximum pivot element, and
exchange rows when it becomes necessary.>> So...my burning
question, in three closely-related parts, are:>> 1) Why is
pivoting necessary? Is it a safeguard against>
divide-by-zero?> 2) If I were to skip pivoting altogether,
why would this> de-stabilize the algorithm?> 3) What
scenarios could come up that would fall apart in the
Fritz>> With regard to part 2>> The choice to pivot or not is
strongly>> in§uenced by the kind of matrix data. When>> using
the Crout method of LU decomposition>> for matrix inversion, I
was able to invert>> some fairly large matrices without
pivoting.>> They were, however, 'lled with random numbers 
in>>
the range [-10., 10.]>> my results were>> size time ave
accuracy>> N millisec>> no pivoting>> 100 60 3. e-15>> 200 440
4. e-14>> 400 6700 7. e-15>> with partial pivoting>> 100 60 1.
e-15>> 200 550 4. e-14>> 400 7080 7. e-15>> accuracy is simply
the norm of input A times>> inverse AI minus unit matrix I.
double>> precision was used on a ordinary Pentium 3 PC>> acc =
norm(A*AI-I)/N*N>> The pivoting made little change, because the
data>>was so nice.>> There are also some realy nasty matrices
like>>the Hilbert where pivoting won't help>> Bill
Shortall>>
===
>>Subject: Re: Q. On Crout's Algorithm For Matrix
student, trying to write a software>> library for handling
matrices that, among other things, 'nds>> determinants
ef'ciently. The best way to do this is by LU>> decomposition,
and the best way to do that (that I know of) is Crout's>>
algorithm.> I stumbled upon Crout's algorithm in _Numerical
Recipes In C: The Art>> of Scienti'c Computing_, which has a
pretty good discussion. I>> understand most of it.>>What is
described in NR (2.3 Performing the LU decomposition) is>>not,
strictly speaking, the Crout's algorithm. NR suggest
the>>following scheme: compute the 1st column of U, then the
1st >>column of L, the 2nd column of U, the 2nd column of U,
etc.>>C.f. G.E.Forsythe, C.B.Moler, Computer Solution of
Linear>>Algebraic Systems, Prentice-Hall, Englewood Cliffs,
N.J., 1967,>>where Crout's algorithm is described using the
other sequence:>>compute 1st row of U, 1st column of L, 2nd
row of U, 2nd column>>of L, etc.>>Provided that the latter was
published prior to the former,>>I think that Crout's method 
is
the row-column one.>>Actually, in that same section 12, Crout
and Dolittle Modi'cations,>>Forsythe&Moler put a note that 
the
algorithm being described is>>in fact the Dolittle method,
while Crout algorithm is a modi'cation>>of the Gauss LU
decomposition, where U is unidiagonal (u(i,i)=1).>>In the
Gauss LU decomposition the main diagonal of L is all 1,>>i.e.
L is unidiagonal.>> Where I get stuck is the concept of the
>pivoting> (Swapping rows,>> when appropriate. This results in
a row-wise permutation of the>> the stability of Crout's
method.> Calculation of the subdiagonal>> elements
(alpha[i][j]) involves dividing by a pivot element. The>>
mathematical discussion of the algorithm just uses the
non-unity>> diagonal elements (beta[j][j]), without regard to
pivoting. I worked>> through it by hand with a random 3x3
matrix, and it seemed to work>> just 'ne. However, the 
example
C code goes out of its way to choose>> maximum pivot element,
and exchange rows when it becomes necessary.> So...my
burning question, in three closely-related parts, are:> 1)
Why is pivoting necessary? Is it a safeguard against>>
divide-by-zero?>> 2) If I were to skip pivoting altogether,
why would this>> de-stabilize the algorithm?>> 3) What
scenarios could come up that would fall apart in the absence>>
Ralston, >A 'rst course in numerical analysis> 2nd ed>
(McGraw-Hill, New York, 1978). It discusses the pivoting
question in> sickening detail.>>I would suggest the basic
works of Jim Wilkinson>>Error Analysis of Direct Methods of
Matrix Inversion,>>J. Assoc. Comput. Mach., 8 (1961)
281-330>>Rounding Errors Algebraic Processes, Prentice-Hall,
Englewood Cliffs,>>N.J., 1963>>The Algebraic Eigenvalue
Problem, Clarendon Press, Oxford, 1965>>I guess there was
nothing new for Gauss-like methods error analysis>>since then;
besides >partial> and >full> pivoting are the terms>>suggested
by Wilkinson himself.>>Roughly speaking (in a non-formal way),
pivoting is a method to>>prevent growth in derived matrices
elements. Why this is important>>(to bound the growth of
matrix elements): a computer solution y>>is an exact solution
of a perturbed linear system (due to a large>>number of
rounoffs, on the order of O(N^3)):>> (A + F)x = b + f>>F is a
perturbation matrix, f is a perturbation in the RHS.>>The norm
of the perturbation matrix F is proportional to>>the growth of
the elements of derived matrices on the forward>>steps of
Gauss method. Therefore, large growth results in a>>solution
of a system that is by far (of by F) deviated from>>the
original system with matrix A, and the answer to the>>question
1. of the OP is yes, but not only a zero division>>safeguard.
The answer to the second question should be more>>or less
obvious as well.>>The idea behind pivoting should be clear
now, and for the>>full theory of the backward error analysis I
am referring to>>the Wilkinson's papers.> Pivoting is
necessary to reduce roundoff error. Any of the iterative>
elimination schemes for solving linear equations or inverting
matrices> involves multiplying a row by a constant and
subtracting it from another> row. >>I think this is inexact:
Gauss-like LU decomposition methods are not iterative,>>on the
contrary they are 'nite methods, that terminate in a 
'nite
number of>>steps even in exact math. Iterative methods in
exact math converge after>>in'nite number of iterations.> It
is the subtraction step that hurts because it frequently
happens> that you subtract numbers that are close in
magnitude. The difference has> many fewer signi'cant 
'gures
of precision than either of the numbers> you were
subtracting. By pivoting and normalizing to the largest
element> in a row one reduces this (unavoidable) loss of
destructive. The actual reason to>>introduce pivoting methods
is to reduce growth in matrix elements,>>not to avoid
D.Goldberg paper >What every Computer Scientist
should>>know...>> form of pivoting involves permuting both
rows and columns, but the resulting> algorithm is far mor
complex than row-wise pivoting alone. Fortunately,> in most
practical cases row-wise pivoting is adequate.>>Not that far
more complex, it would be rather more slow, since it
involves>>O(n^3) comparisons, not O(n^2) as in the partial
pivoting case.>>The advantage of the full pivoting over the
partial pivoting is that>>it ensures a better upper bound for
the growth in the matrix elements>>in the elimination process,
and therefore for the smaller norm of the>>perturbation
matrix.>>Let A(k) be the major submatrix of order k on k-th
forward step>>of the elimination process (in rows and columns
k through n) and>>a(k) = max_{i,j>k}|A_{i,j}(k)|, and let the
growth factor be >>g_n(A) = max_{0<=k<=n-1}a(k)/a(0). It has
been shown that the upper>>bound for g_n(A) in the partial
pivoting scheme is g_n(A) < 2^(n-1).>>It is large and there
are matrices for which this bound can be >>actually reached,
e.g. (a well known example!):>> 1 0 0 0 . . 0 1>>-1 1 0 0 . .
0 1>>-1 -1 1 0 . . 0 1>>-1 -1 -1 1 . . 0 1>> . . . . . . . .>>
. . . . . . . .>>-1 -1 . . . -1 1 1>>-1 -1 . . . . -1 1>>The
fact is that without a pivoting there is *no* upper bound for
the>>growth factor, it can be arbitrary large. Therefore the
answer the>>question 3. is there is a possiblity to obtain an
arbitrary astray>>solution in the absence of a pivoting.>>For
the full pivoting the upper bound is smaller and
asymptotically>>is >>g_n(A) <= f(n) prop Cn^(1/4log(n))>>In
particular, for the above matrix the full pivoting gives
g_n(A) <= 2>>for any n.>>Rgds,>>Andrew>> -- > Julian V.
Noble> Professor Emeritus of Physics>
http://galileo.phys.virginia.edu/~jvn/>>God is not willing
to do everything, and thereby take away our free wiil> and
that share of glory that rightfully belongs to
ourselves.>> -- N. Machiavelli, >The
Prince>.>>
===
>>Subject: Does 0 have a signi'cant digit?>>All
you.>>
===
>>Subject: Re: algorithms for skeleton tracing and
loop extraction> Do you know of any >standard> algorithms
for skeleton tracing and loop> extraction? Here are the
details:>> What I have:> I have a binary image that
contains a skeleton of a contour. This> skeleton can be
thought of as a collection of links connecting feature>
points of three types: end points, T- and X-intersections.
This> skeleton may contain loops of> various degrees of
complexity, e.g. a simple loop where a link starts> and ends
at the same FP or a complex loop where several FPs are>
connected into a single loop.>> What I need:> I need to
trace this contour, i.e. to turn a 2d image into a sequence>
of skeleton pixels going from one FP to another in a certain
order. I> also need to detect loops. This implies that when
traced, all links> belonging to the same loop should be
listed successively without any>gaps>.>>I think you're
asking for two separate things here. The 'rst>>is a way of
turning your image into a collection of links>>between FPs.
The second is the >tracing and loop
extraction>.>>Unfortunately, I'm not sure exactly what you
want the second>>part to do. What it sounds like you're 
asking
for isn't always>>possible; some patterns of linkage 
*can't* be
traced out in>>a single sequence of adjacent pixels, and even
when they can>>it isn't always possible to keep all links 
that
belong to the>>same loop together. I'll come back to that
later; 'rst of all,>>here's how to turn pixels 
into links. I'm
sure there are standard>>algorithms for this. But, since I
happen not to know them :-),>>here's a suggestion.>>0. Note
which >set> pixels are adjacent. (The time it takes>> to do
this is linear in the number of >set> pixels.)>>1. De'ne a
>pre-FP> to be any >set> pixel that has at least>> two
adjacent >set> neighbours.>>2. Coalesce adjacent pre-FPs. (The
time it takes to do this>> is more or less linear in the number
of pre-FPs.) Call>> the equivalence classes >vertices>.>> It
may produce better results if you only coalesce pre-FPs>> when
they *overlap*.>>3. De'ne an >edge pixel> to be any set pixel
that isn't>> a pre-FP. Remark: an edge pixel can have at 
most
two>> neighbours.>>4. Coalesce adjacent edge pixels. (The time
it takes to do>> this is more or less linear in the number of
edge pixels.)>> Call the equivalence classes >edges>.>>
Theorem: an edge has either 2 >ends>, or none; in the latter>>
case it forms a simple closed loop. It's possible, and not>>
too painful, to do the coalescing in such a way that you>>
identify the ends (if any) and put the pixels into a linear>>
sequence joining the ends (or a circular sequence including>>
all of them, if there are no ends).>>5. Now look for
adjacencies between vertices and edges.>> Note that it is
impossible for two vertices to be>> adjacent, or for two edges
to be adjacent. (But if>> you adopt the more conservative
pre-FP coalescing>> policy mentioned above, then two vertices
can be adjacent;>> in that case, invent a new edge joining
them.)>>6. Choose a single pixel belonging to each vertex,
and>> 'nd a minimum-length path joining it to each other>>
pixel in that vertex. (The time it takes to do this>> is
linear in the number of vertices, provided the>> number of
pixels in a vertex is bounded. The constant>> of
proportionality is small provided the number of>> pixels in a
vertex is small. Which I bet it is.)>>OK. Now you've got
enough information to do your loop>>extraction and to list the
pixels in each loop (or,>>more generally, in any path through
the vertices and>>edges) in a sensible order. (Start at the
chosen pixel>>in the starting vertex. Follow the shortest path
from>>there to the nearest point in the next edge. Then
follow>>the edge's pixels in order. Then follow the shortest
path>>from there to the chosen pixel in the next vertex.
And>>so on.) What about tracing and loop extraction?>>I said,
earlier, that it isn't always possible to do>>what you seem 
to
be asking for. To see why, suppose>>there are four FPs A,B,C,D,
with links between all>>pairs. (To draw this, put A,B,C in a
triangle, D in>>the middle.) Then:>> - There is no way to draw
all the links in a single>> traversal without breaks.>> - Given
any pair of links, there's a loop that includes>> both. But 
you
can't list all the links so that there>> are no >gaps> in 
any
of these loops.>>I suspect I've just misunderstood what 
you're
asking for.>>Can you be a bit more speci'c?>>-- >>Gareth
McCaughan>>.sig under construc>>
===
>>Subject: Mutli-parameter
border search>>I have a simulation that takes in quite a few
parameters (say around>>10). The data output can be simpli'ed
to just one point on a>>two-dimensional chart. >>I am
interested in is the border of that plot (i.e. the
outlying>>points), not the points in the middle. The
relationship between the>>input parameters and the output
points is not straightforward, though>>not random. What is the
best way to >move> around the border of the>>plot when running
the simulation (This will save much simulation time,>>as I'd
avoid running most simulations that>>generate points around
the center of the plot)?>>What's the best approach? Is it
using GAs - or a dimensional mapping>>method such as
Neldor-Mead? I'm a bit lost with this, so any>>suggestions
===
>>Subject: L1
linear 't>>Please point me to an ef'cient 
routine for linear
L1 't, that is,>>'nding a minimizer of x |-> 
norm(b-A*x, 1)
with a given column vector b>>and a (generally speaking)
rectangular matrix A. Is there a routine about>>as
well-researched as least-squares routines?>>I am a MATLAB
user, so a routine in this system would be most welcome,
but>>a description in mathematical formulas would be
satisfactory.>>I tried Google, but perhaps didn't set the 
key