Cholesky solution of matrices
Post by Himself on Apr 13th, 2017, 10:13am
Code: Cholesky App
:
Andre-Louis Cholesky (1885 - 1918) : French military officer and mathematician.
Test data extracted from J. C. Nash 1979 publication: who used a HP 9830 workstation.
:
The Cholesky algorithm is here tittivated for precise Engineering Surveying applications,
and others busying themselves with multi-parameter least squares calculations.
:
The advent of the *FLOAT 80 function enhances the validity of the "most likely" solution
for the sought value of the parameters ( to 19 significant figures no less ).
:
Modern surveying instruments called Total Stations can provide sub-second and sub-millimetre values.
This requires computational methods to match ongoing technological developments and engineering challenges.
:
Run under Windows 10 Pro (1607): BB4W v 6.10a : Computer HP Z400.
:
Himself
Re: Cholesky solution of matrices
Post by Richard Russell on Apr 14th, 2017, 09:11am
on Apr 13th, 2017, 10:13am, Himself wrote:| The Cholesky algorithm is here tittivated for precise Engineering Surveying applications, and others busying themselves with multi-parameter least squares calculations. |
|
The code box in your post is empty!
I don't know the Cholesky algorithm, but BBC BASIC is exceptionally well suited to least-squares-fitting calculations because of its built-in matrix arithmetic, especially the dot-product operator. One of the example programs supplied with BBCSDL, polyfit.bbc, illustrates this technique.
Richard.
Re: Cholesky solution of matrices
Post by Himself on Jun 11th, 2017, 1:56pm
:
REM Cholesky App
:
REM Andre-Louis Cholesky (1885 - 1918) : French military officer and mathematician.
REM Test data extracted from J. C. Nash 1979 publication: who used a HP 9830 workstation.
REM
REM The Cholesky algorithm is here tittivated for precise Engineering Surveying applications,
REM and others busying themselves with multi-parameter least squares calculations. It is
REM space-saving and speedy.
REM
REM The advent of the *FLOAT 80 function enhances the validity of the "most likely" solution
REM for the value of the parameters ( to 19 significant figures no less ).
REM
REM Modern surveying instruments called Total Stations can provide sub-second and sub-millimetre measurement values.
REM This requires computational methods to match ongoing technological developments and engineering challenges,
REM without the arithmetical burden, and limitations, of the matrix inversion approach to the sought-after solution.
REM
REM Run under Windows 10 Pro (1607): BB4W v 6.10a : Computer HP Z400, Quad core 2.40 Xeon.
REM
:
*FLOAT 80
*FONT Courier,16,B
:
INSTALL @lib$+"FNUSING"
:
(Start)
INPUT " For this demonstrationn enter 5 as the number of unknown parameters for this test";n%
IF n% <> 5 THEN
SOUND 1,-15,52,-1
SOUND 1,-15,56,50
WAIT 50
SOUND &11,-15,100,10
CLS
GOTO (Start)
ENDIF
:
: REM Normal equations n% x n% symmetric, positive definite.
:
DIM a(n%*(n%+1)/2) : REM Vector storage (n*(n+1))/2 where n equals the
DIM x(n%) : REM order of the matrix.
DIM q%(n%)
: REM Vector storage of the lower triangle of the normal equations.
a(1)=18926823
a(2)=6359705 : a(3)=2164379
a(4)=10985647 : a(5)=3734131 : a(6)=6445437
a(7)=3344971 : a(8)=1166559 : a(9)=2008683 : a(10)=659226
a(11)=14709 : a(12)=5147 : a(13)=8859 : a(14)=2926 : a(15)=13
:
: REM Vector storage of equations' right-hand side.
:
x(1)=5937938 : x(2)=2046485 : x(3)=3526413 : x(4)=1130177 : x(5)=5003
:
FOR j% = 1 TO n% : REM Computing the decomposition of the normal equations.
q%=((j%*(j%+1)) DIV 2)
IF j%>1 THEN
FOR i% = j% TO n%
m% = ((i%*(i%-1)) DIV 2)+j% : s=a(m%)
FOR k% = 1 TO (j%-1) : s=s-a(m%-k%)*a(q%-k%) : NEXT k%
a(m%)=s
NEXT i%
ENDIF
: REM Cholesky factor.
s=SQR(a(q%))
FOR i% = j% TO n%
m%=((i%*(i%-1)) DIV 2)+j%
a(m%)=a(m%)/s
NEXT i%
PRINT
NEXT j%
:
VDU 30 : CLS : PRINT SPC(5);"Cholesky factor output vector, in row order."
FOR i% = 1 TO n%
FOR j% = 1 TO i%
PRINT SPC(5);a(j%)
NEXT j%
PRINT
NEXT i%
: REM Forward substitution commences.
x(1) = x(1)/a(1)
q%=1
FOR i%=2 TO n%
FOR j%=1 TO (i%-1)
q%=q%+1 : x(i%)=x(i%)-a(q%)*x(j%)
NEXT j%
q%=q%+1x(i%)=x(i%)/a(q%)
NEXT i%
: REM Backward substitution from here.
x(n%)=x(n%)/a(n%*(n%+1)/2)
FOR i%=n% TO 2 STEP -1
q%=i%*(i%-1)/2
FOR j%=1 TO (i%-1)
x(j%)=x(j%)-x(i%)*a(q%+j%)
NEXT j%
x(i%-1)=x(i%-1)/a(q%)
NEXT i%
:
PRINT SPC(5);"Determinate solution for the 5 unknowns in aproximately 30 milliseconds." : PRINT
FOR i% = 1 TO n%
PRINT SPC(5);"x(";STR$(i%);") = ";FNusing("###.##################",x(i%))
NEXT i%
:
:
END!!
:
:
Re: Cholesky solution of matrices
Post by Richard Russell on Jun 11th, 2017, 3:49pm
on Jun 11th, 2017, 1:56pm, Himself wrote:| REM The advent of the *FLOAT 80 function enhances the validity of the "most likely" solution |
|
Just to note that *FLOAT 80 has no effect at all in this program, since it affects only floating-point indirection and data files, and the program uses neither.
Richard.