# GPS: minimizing sum of squares algorithm

### [1/14] from: jjmmes:yaho:o:es at: 8-Nov-2002 19:41

Has anybody written this algorithm in Rebol ?
I need this to georeference gps data, i.e. find a, b,
c, d, e, f below given three calibration points
x = a + b lat + c lon
y = d + e lat + f lon
It's available as a math library in most languages but
maybe some rebol guru has it too !
Regards

### [2/14] from: jan:skibinski:sympatico:ca at: 8-Nov-2002 14:22

Jose,
I could certainly help you with this but I would
need to know more about the shape of the functions b c e f
in order to model them appriopriately. For example,
depending on what those functions represent a polynomial
approximation might be of no use or it might be
quite good actually.
A reference to existing georeference algorithms could
be a good starting point.
Regards,
Jan
jose wrote:

### [3/14] from: tomc:darkwing:uoregon at: 8-Nov-2002 12:12

Hi Jose
I implemented some linear alrerbra functions a few years back.
I will see if I can dig them up over the weekend
On Fri, 8 Nov 2002, [iso-8859-1] jose wrote:

### [4/14] from: jjmmes:yaho:o:es at: 9-Nov-2002 16:23

a,b,c,d,e,f are constants.
I basically want to solve a set of 6 equations with
six unknowns, using one of the existing linear algebra
algorithms. The algorithm implementation I'm looking
for is also known as RMS, root mean squares but
couldn't find a good algorithm description in the web
Regards
Jose
--- Jan Skibinski <

**[jan--skibinski--sympatico--ca]**> escribió: > Jose,### [5/14] from: jan:skibinski:sympatico:ca at: 9-Nov-2002 11:01

Jose,

> a,b,c,d,e,f are constants.
> I basically want to solve a set of 6 equations with

<<quoted lines omitted: 4>>

> Regards
> Jose
Ah, then this is quite simple. I thought you used Rebol syntax
with equal signs put there by mistake instead of ( :), while
in fact you wrote mathematical equations with implied *
between 'b and 'lat, etc. Sorry for that.
I'll post a solution later today.
Jan

### [6/14] from: jan:skibinski:sympatico:ca at: 10-Nov-2002 3:40

Jose,
As promised, here is a solution to your GPS minimization
problem, or rather a part of it. The only function I am
attaching here is the main GPS function. I have not tested
it, so it might be buggy, therefore threat it as the outline
only. I developed and tested today various pieces that
the gps function is calling, including a solver for a system
of linear equations. If you want to use them just ask for
it and I will send it to you. I did not plan to develop any
linear algebra package at this time and whatever I developed
today is not optimized and quite rough around the edges.
For this reason I am not publicly posting it to the library.
But you can use any other solver, or any other package with
the gps function below.
All the best,
Jan
comment {
Georeference GPS data modeling:
------------------------------
Find six constants relating [x y] coordinates
and [long lat] coordinates, which we abbreviate
here as [u v]:
x = a1 + a2 u + a3 v
y = b1 + b2 u + b3 v
To find a1, a2, a3, b1, b2, b3 we need at least 3 sets
of data for x, y, u and v - the more the merrier:
[x1 x2 x3 ... xn]
[y1 y2 y3 ... yn]
[u1 u2 u3 ... un]
[v1 v2 v3 ... vn]
Minimization of the least square error F:
F = Sum (xi - a1 - a2 ui - a3 vi)**2
+ (yi - b1 - b2 ui - b3 vi) ** 2
Six stationary conditions for extremum of F:
dF/da1 = 0, dF/da2=0, dF/a3=0,
dF/db1=0, dF/db2=0, dF/db3=0
Simplification:
-------------
I assumed below the same average variance for all the
samples,which is good enough if your data is reasonably
correct, otherwise the model would need to deal with
individual variances, and the description would not be
as simple as this one.
And because the variance is constant it will not appear
anywhere in the equations that follow.
Since variables x and y are independent the system of
six equation separates itself into two system of three
equation each.
We have to separately solve two matrix equations,
M a = p
M b = q
The matrix M is:
| n | Sum ui | Sum vi |
| | | |
| Sum ui | Sum ui ui | Sum ui vi |
| | | |
| Sum vi | Sum ui vi | Sum vi vi |
Vector p-transposed = [Sum xi, Sum xi ui, Sum xi vi]
Vector q-transposed = [Sum yi, Sum yi ui, Sum yi vi]
Internally, we will represent matrix M as block of
three block-columns.
}
gps: func [
{A block of two blocks [[a1 a2 a3][b1 b2 b3]]
containing estimated model parameters, computed
from 4 sets of input data 'xs 'ys 'us 'vs
}
xs [block!] {[x1 x2 .... xn]}
ys [block!] {[y1 y2 ..... yn]}
us [block!] {[u1 u2 ..... un]}
vs [block!] {[v1 v2 ..... vn]}
/local n m p q result [block!]
][
simple-sum: func[x][foldl :+ 0 x]
n: length? xs
m: copy []
append/only m reduce [n simple-sum us simple-sum vs]
append/only m reduce [simple-sum us product-vv us us product-vv
us vs]
append/only m reduce [simple-sum vs product-vv us vs product-vv
vs vs]
p: reduce [simple-sum xs product-vv xs us product-vv xs vs]
q: reduce [simple-sum ys product-vv ys us product-vv ys vs]
result: copy []
append/only result solve m p
append/only result solve m q
result
]
comment {
'product-vv stands for scalar inner product of two vectors.
'Solve uses my home grown algorithm based on Gram-Schmidt
orthogonalization procedure. A typical approach here would be
to use the LU decomposition technique with pivoting.
}

### [7/14] from: jjmmes::yahoo::es at: 10-Nov-2002 18:11

Jan,
Many thanks for your help. Can you pls send me your
functions so that I can test your gps functions
I find that there are very few math related scripts
available in rebol, so your scripts are most valuable!
Thanks
--- Jan Skibinski <

**[jan--skibinski--sympatico--ca]**> escribió: > Jose,### [8/14] from: jan:skibinski:sympatico:ca at: 10-Nov-2002 15:14

Jose and all,

> Many thanks for your help. Can you pls send me your
> functions so that I can test your gps functions
>
> I find that there are very few math related scripts
> available in rebol, so your scripts are most valuable!
>

Jose, I just sent you the script you requested. Here are
just few general comments about linear algebra scripts.
I am sure several people here developed their
own Linear Algebra packages in Rebol.
Tom Conlin for one. I bet Joel Neely has one too.
I have developed several Haskell modules that deal
with some aspects of Linear Algebra, and go way beyond
the basic stuff, such as LU decomposition and a solver
of systems of linear equations. I needed support for
nonsymmetric complex matrices, eigenvalues, eigenvectors,
and all that stuff. I have also developed a Quantum Vector
module which supports the <bra| |ket> Dirac's notation,
and several application modules from Quantum Mechanics,
one of which is a prototype of
"Haskell Simulator of Quantum Computer".
I did no plan however to port this stuff to Rebol, knowing
well that this is not in a general interest of Rebolites.
But your GPS problem forced me to port some of
the stuff to Rebol.
Here are few of my observations about my yesterday's
porting process:
------------------------------------------------------
The process was rather easy, because of the clarity
of the Haskell syntax. I did not spend much time on optimization
and I used the same concepts in Rebol as in Haskell.
So, with few exceptions, the recursive stuff in Haskell,
remains recursive in Rebol. The original functions were
by defnition "indexless" - they came from my web
article "Indexless Linear Algebra Algorithms."
The resulting translation is almost indexless too.
The list comprehension notation , such as:
p = [bra_ket u ak | ak <- a]
were easy to translate to
p: copy []
foreach ak a [insert tail result (bra-ket u :ak)]
However, since Haskell is lazy the expressions
are not evaluated until they are needed. So if during
recursion, some items such as 'u cannot be evaluated
because say u: first a but 'a is already empty, nothing
wrong happens in Haskell because of its laziness.
Not so in Rebol.
The solution is to use functions, even if they represent
constants. So the above 'p would be wrapped
inside a function instead.
p: func [
/local result
][
result: copy []
foreach ak a [insert tail result (bra-ket u :ak)]
result
]
Haskell has several standard words that are identical
as in rebol but with completely different
semantics: 'head means 'first, 'tail means 'next, etc.
But one can easily get used to such things.
Few things have bitten me though. I am still not used
to certain "in-place" behaviour of certain functions, such
as 'replace, or worst of all 'reverse, which returns empty
block. This bites badly in pipelines of expressions.
Regards,
Jan

### [9/14] from: jjmmes::yahoo:es at: 10-Nov-2002 21:36

Can you pls send a link to your web article
Indexless Linear Algebra Algorithms
?
The link at numeric-quest seems broken !
Thanks again
--- Jan Skibinski <

**[jan--skibinski--sympatico--ca]**> escribió: >### [10/14] from: tomc:darkwing:uoregon at: 10-Nov-2002 13:55

Hi Jose,
Hmmm
I recall "root mean square" from Calculus not Linear
and it was the sqrt of the area under f(x^2) between a couple of points
if this is the correct interpretation I am not seeing how it will help you
find your constants.
I think you are closer to your mark in the subject line
and may be looking for a "linear-least-squares" fit.
(or if you are less lucky a "non-linear-least-squares" fit)
disclaimer: I have only been through standard undergraduate classes and
have not used this stuff in several years
(winter 99, and never in the real world).
if I were you I would try to view what I needed
as 2 problems each with 3 unknowns since your
x and y equations seem independent and the amount of work will
normally square with the number of unknowns,
that is 2*3^2 < 6^2
looking at just the x half,
a few gps samples you could have in something similar to your notation
1*a + lat1*b + lon1*c = x1
1*a + lat2*b + lon2*c = x2
1*a + lat3*b + lon3*c = x3
or in matrix form
| 1 lat1 lon1 | |a| |x1|
| 1 lat2 lon2 | * |b| = |x2|
| 1 lat3 lon3 | |c| |x3|
in a more standard notation this is typically represented as
A x = b
the linear-least-squares solution 'x' to Ax = b is
x = (A_t * A)^-1 * A_t * b
where A_t is the transpose of A
(transpose means the indices for a row and column are swapped
or that a matrix is reflected about its diagonal)
and the (A_t A)^-1 means the inverse of the matrix that results from
multiplying the transpose of A by A (order matters)
more specifically it is the matrix you would have to multiply (A_t * A)
by to get the "Identity matrix" where the identity matrix is
a square matrix with 1's from the top-left to lower-right and 0's
elsewhere
a matrix must be square to have an inverse
There is no guarantee a square matrix has an inverse.
I am not going into inverting a matrix because there are alternatives.
all in all it is pretty expensive to go this way
but at least it is straight forward if you want to plow through it.
If you are going to do allot of these there is a shortcut which
stems from the rearrangement
A_t A x = A_t b
where the two sides of the equation can be expressed as
Augmented Matrices
going back to something that resembles your notation translated to matrix
form we get
|| 1 1 1| |1 lat1 lon1| : a| | 1 1 1 : x1|
||lat1 lat2 lat3| * |1 lat2 lon2| : b| = |lat1 lat2 lat3 : x2|
||lon1 lon2 lon3| |1 lat3 lon3| : c| |lon1 lon2 lon3 : x3|
where the [a b c] & [x1 x2 x3] vectors are augmenting the square matrices
these augmented matrices (representing your simultaneous equations)
may be solved by adding multiples of the rows in a matrix to one another.
to end up with 1s, 0s, and solutions (if there are any)
there are methods by Gauss, and an improvement called "Gauss-Jordan"
which have been around for hundreds of years.
in the end you end up with identity matrices augmented with the solutions
|1 0 0 : i*a| |1 0 0 : p|
|0 1 0 : j*b| = |0 1 0 : q|
|0 0 1 : k*c] |0 0 1 : r]
so a = p/i
b = q/j
c = r/k
there is a function "naive-gauss" in code I wrote to double check
the longhand results required in my homework.
(the professor would not even allow a four function calc in the class)
and was not finished beyond the point where it served my immediate
purpose
http://darkwing.uoregon.edu/~tomc/core/math/matrix.r
that should work to reduce your augmented equations.
there is another "gauss-jordan" which would be better if it did not have
a bug I never got around to tracking down
note:
the "comp" function in the file was derived from Eric Long's compare
best of luck
On Sat, 9 Nov 2002, [iso-8859-1] jose wrote:

### [11/14] from: jan:skibinski:sympatico:ca at: 10-Nov-2002 16:48

Hi Jose,

> Can you pls send a link to your web article
> "Indexless Linear Algebra Algorithms" ?
>
> The link at numeric-quest seems broken !
>
> Thanks again
>

I know, I broke it. The site is still defunct. But you can retrieve
its past contents from the Wayback Machine. Go to
www.archive.org, type in a name of a site you are interested
in, and voila! Isn't that neat?
While in Numeric Quest, select page with Haskell modules.
The recommended way is to openit in a separtate window,
so you will break out of the stupid frames.
All so-called "literate modules", which look like HTML
pages, used to be directly executable in Hugs interpreter.
Just rename a file giving it a literate extension "lhs" rather
than html and you should be able to run it - provided that
Wayback Machine did not do too much damage while
reformatting the pages - like substituting ">" by >,
or by adding java scripts to them. But the Orthogonals.html
(from the last archive) seems to be accepted by Hugs
in spite of all that reformatting.
Enjoy,
Jan

### [12/14] from: g:santilli:tiscalinet:it at: 11-Nov-2002 11:31

Hi Jan,
On Sunday, November 10, 2002, 9:14:54 PM, you wrote:
JS> Few things have bitten me though. I am still not used
JS> to certain "in-place" behaviour of certain functions, such
JS> as 'replace, or worst of all 'reverse, which returns empty
JS> block. This bites badly in pipelines of expressions.
I don't know why REVERSE returns the tail of the block. You can
use HEAD REVERSE in expressions, or patch it like this:
use [_reverse] [
_reverse: :reverse
reverse: func [
"Reverses a series."
value [series! tuple! pair!]
/part "Limits to a given length or position."
range [integer! series!]
] [
head either part [_reverse/part value range] [_reverse value]
]
]
Regards,
Gabriele.
--
Gabriele Santilli <

**[g--santilli--tiscalinet--it]**> -- REBOL Programmer Amigan -- AGI L'Aquila -- REB: http://web.tiscali.it/rebol/index.r### [13/14] from: jan:skibinski:sympatico:ca at: 11-Nov-2002 11:33

Hi Gabriele,

> I don't know why REVERSE returns the tail of the block. You can
> use HEAD REVERSE in expressions, or patch it like this:

<<quoted lines omitted: 9>>

> ]
> ]
Thanks, that's interesting application of 'use! But I'd rather stay
away
from patching and just get used to the standard behaviour of 'reverse.
Jan

### [14/14] from: jjmmes:ya:hoo:es at: 11-Nov-2002 21:21

Hi Tom,
Thanks for helping me out. I think the general
approach to solving this is with Jan's code, specially
when you want to acount for 4+ data points.
It turns out that with 3 data points I only have 2
sets of independent equations, u: f(x,y) and v:
f'(x,y), so I only have to solve one set of 3
equations with 3 unknowns and can get the unknows
expressed as a function of u,v .
Regards,
Jose
--- Tom Conlin <

**[tomc--darkwing--uoregon--edu]**> escribió:> Hi Jose,
> Hmmm

<<quoted lines omitted: 153>>

> > > > c, d, e, f below given three calibration
> points
=== message truncated ===

Notes

- Quoted lines have been omitted from some messages.

View the message alone to see the lines that have been omitted