Fortran - cholesky decompostion and inverses of upper diagonal matrices

Author

cholesky decompostion and inverses of upper diagonal matrices

student

2005-06-07, 4:02 pm

i am currently using the routine available with numerical recipes to do
cholesky decomposition? i would like to get inputs from people about
the computation speed of that scheme.

also, are there any fast routines designed specifically to invert upper
diagonal matrices? these upper diagonal matrices are obtained by
cholesky decomposition of symmetric positive definite matrices.

jon

2005-06-07, 8:58 pm

I don't know about speeds, but I believe the code provided by Alan
Miller at

http://lib.stat.cmu.edu/apstat/274
or
http://users.bigpond.net.au/amiller/
(look at code for lsq.f90)

includes a routine "inv" that does just that --> (invert upper
diagonal matrices? these upper diagonal matrices are obtained by
cholesky decomposition of symmetric positive definite matrices. )

You can tweak those algorithms as you wish, but there may be faster
algorithms to obtain a solution.

HTH

Salvatore

2005-06-08, 8:58 am

student wrote:
> i am currently using the routine available with numerical recipes to do
> cholesky decomposition? i would like to get inputs from people about
> the computation speed of that scheme.
>
> also, are there any fast routines designed specifically to invert upper
> diagonal matrices? these upper diagonal matrices are obtained by
> cholesky decomposition of symmetric positive definite matrices.

Use a combination of LAPACK with the ATLAS BLAS, from
http://www.netlib.org/
the routines you need are DPOTRF and DPOTRI. And, you mean triangular,
right? not diagonal.

Assuming, of course, that you REALLY need an explicit inverse,
something that should always be regarded with suspicion, 99 out of 100
times it is simply not true.

Hope this helps
Salvatore Filippone

David Frank

2005-06-08, 4:00 pm

"Salvatore" <sfilippone@uniroma2.it> wrote in message
news:1118227117.882976.229080@o13g2000cwo.googlegroups.com...

<snip>
>
> Assuming, of course, that you REALLY need an explicit inverse,
> something that should always be regarded with suspicion, 99 out of 100
> times it is simply not true.
>
> Hope this helps
> Salvatore Filippone
>

How about you or someone proving lapack solution can beat my code that uses
inverse
that was just posted in "wanted fast routines" topic

http://home.earthlink.net/~dave_gemini/linear.f90

http://home.earthlink.net/~dave_gemini/linear.exe for comparison between
windows solutions

Victor Eijkhout

2005-06-08, 4:00 pm

David Frank <dave_frank@hotmail.com> wrote:

> How about you or someone proving lapack solution can beat my code that uses
> inverse

If all you need is a system solution, then computing the inverse wastes
operations and is numerically less stable.

Maybe your code runs faster than mine because you're a better coder. If
you'd coded the right algorithm yours would run even faster!

V.
--
email: lastname at cs utk edu
homepage: www cs utk edu tilde lastname

Rich Townsend

2005-06-08, 4:00 pm

David Frank wrote:
> "Salvatore" <sfilippone@uniroma2.it> wrote in message
> news:1118227117.882976.229080@o13g2000cwo.googlegroups.com...
>
> <snip>
>
>
>
> How about you or someone proving lapack solution can beat my code that uses
> inverse
> that was just posted in "wanted fast routines" topic

I just did -- see the "wanted" thread.

:)

cheers,

Rich

David Frank

2005-06-08, 4:00 pm

"Rich Townsend" <rhdt@barVOIDtol.udel.edu> wrote in message
news:d87arl$t2o$1@scrotar.nss.udel.edu...
> David Frank wrote:
>
> I just did -- see the "wanted" thread.
>
> :)
>
> cheers,
>
> Rich

Hi Rich,
try my latest and greatest, which gives better results for coef =50 and
outputs absolute error..

http://home.earthlink.net/~dave_gemini/linear.f90

Rich Townsend

2005-06-08, 8:58 pm

David Frank wrote:
> "Rich Townsend" <rhdt@barVOIDtol.udel.edu> wrote in message
> news:d87arl$t2o$1@scrotar.nss.udel.edu...
>
>
>
> Hi Rich,
> try my latest and greatest, which gives better results for coef =50 and
> outputs absolute error..
>
> http://home.earthlink.net/~dave_gemini/linear.f90
>
>
>

OK:

linear : 0.5407 +/- 3.8e-4
linear_lapack : 0.5364 +/- 3.8e-4

Note that the standard deviations are rather larger than in my previous
posts -- I was making a mistake calculating them. But the difference
between the two is still large at 10 sigma.

How about you write a test code where the timing relates to the matrix
inversion only?

cheers,

Rich

Janne Blomqvist

2005-06-08, 8:58 pm

In article <d87dj2$h1$1@scrotar.nss.udel.edu>, Rich Townsend wrote:
> OK:
>
> linear : 0.5407 +/- 3.8e-4
> linear_lapack : 0.5364 +/- 3.8e-4
>
> Note that the standard deviations are rather larger than in my previous
> posts -- I was making a mistake calculating them. But the difference
> between the two is still large at 10 sigma.
>
> How about you write a test code where the timing relates to the matrix
> inversion only?

Here goes. I split up Davids benchmark, and added another one, as well
as converting it to f95 so that I could compile it with Lahey. Sorry
for the large number of modules, I basically cut-and-glued together
stuff that was lying around on my disk without regard for making a
small self-contained module.

Well, since what separates the wheat from the chaff in linear solvers
really is the accuracy (as someone said, what good is computing the
wrong answer really fast), I added a test for a somewhat more
ill-conditioned matrix than Davids matrix, namely the well-known
Hilbert matrix. For a matrix sufficiently big so that the timings can
actually be measured:

% ./testlinsolve
Hilbert matrix of size 500 is ill-conditioned.
rcond: 4.740663050630338E-21
Condition number: 2.109409568492821E+20
la_gesv solved system in 0.1499999999999999 seconds.
Average cumulative error is: 1577.575979219314
David Frank solved system in 2.570000000000000 seconds.
Average cumulative error is: 4.358666707056873E+17

We see that la_gesv (which uses LU) is more than an order of magnitude
faster. Well, even LU can't get reasonable result with a Hilbert
matrix this size, but for size 11 (size 12 happens to be about the
maximum la_gesv can handle in this case) we see:

% ./testlinsolve
Hilbert matrix of size 11 is well-conditioned.
rcond: 1.913805532331239E-15
Condition number: 522519129089298.4
la_gesv solved system in 0.000000000000000E+00 seconds.
Average cumulative error is: 3.605202376467460E-03
David Frank solved system in 0.000000000000000E+00 seconds.
Average cumulative error is: 18642906443420.82

While the runtimes are neglible in both cases, Franks code has some
serious issues with accuracy.

Here is the code (in addition you need blas, lapack, lapack95):

module conf

implicit none

! At least IEEE 754 double precision.
integer, parameter :: wp = selected_real_kind(15, 307)

end module conf

! ----------------------
module dfrank_linear
use conf

implicit none

contains

! dat(0 = dependent variable
! dat(1:NC = independent variables
subroutine run_frank_test ()
integer,parameter :: NS=100000 ! #data sets
integer,parameter :: NC=50 ! #coef for independent var.
real(wp) :: dat(0:NC,NS), coef(0:NC)
integer :: i
real(wp) :: time1,time2, test(0:NC) = (/ (i, i=nc+1,1,-1) /)
integer :: n
real(wp) :: a(nc+1, nc+1), b(nc+1)

test(1:nc) = test(1:nc)/100
call random_seed()
call random_number(dat)
do n = 1,NS ! embed test coef within random data sets
dat(1:NC,n) = dat(1:NC,n)-0.5d0 + test(1:NC)
dat(0,n) = test(0) + sum(dat(1:NC,n)*test(1:NC))
end do
call make_frank_matrix (dat, a, b)
call cpu_time(time1) ! start benchmark time
coef = linear(a, b, NC)
call cpu_time(time2) ! stop benchmark time
write (*,91) 'e.g. pentium 2.8ghz, CVF, time = ',time2-time1
write (*,91) 'coef recovered = ',coef(0:4), coef(nc)
write (*,91) 'coef avg error = ',SUM(abs(coef-test))/(NC+1)
91 format (a,99(f0.4,1X))
end subroutine run_frank_test

subroutine make_frank_matrix (dat, a, b)
real(wp), intent(in) :: dat(0:,:)
real(wp), intent(out) :: a(:, :), b(:)
real(wp) :: x
integer :: i, j, k, n, nc, ns

nc = size (dat, 1) - 1
print *, nc
ns = size (dat, 2)
print *, ns

a = 0. ; b = 0.
a(1,1) = ns
do n = 1,ns
do k = 1,nc
x = dat(k,n)
j = k+1
a(j,1) = a(j,1) + x
a(1,j) = a(j,1)
do i = 2,nc+1
a(i,j) = a(i,j) + x * dat(i-1,n)
end do
end do
b(1) = b(1) + dat(0,n)
do k = 1,nc
b(k+1) = b(k+1) + dat(k,n) * dat(0,n)
end do
end do
end subroutine make_frank_matrix

! -----------------------
function linear(a, b, nc) result (coef)
implicit none
integer :: nc
real(wp) :: coef(0:nc)
real(wp) :: a(nc+1,nc+1), b(nc+1), temp(nc+1), ai(nc+1,nc+1)
real(wp) :: c, d
integer :: i, j, k, m, ipvt(nc+1)

! - - - invert 'a' matrix - - -
ipvt = (/ (i, i=1,nc+1) /)
DO k = 1,nc+1
m = k-1+MAXLOC(ABS(a(k:nc+1,k)),dim=1)
IF (m /= k) THEN
ipvt( (/m,k/) ) = ipvt( (/k,m/) )
a( (/m,k/),:) = a( (/k,m/),:)
END IF
d = 1/a(k,k)
temp = a(:,k)
DO j = 1,nc+1
c = a(k,j)*d
a(:,j) = a(:,j)-temp*c
a(k,j) = c
END DO
a(:,k) = temp*(-d)
a(k,k) = d
END DO
ai = a(:,ipvt)
! - - - use inverse to evaluate coef - - -
do k = 0,nc
coef(k) = ai(k+1,1) * b(1)
do i = 1,nc
coef(k) = coef(k) +ai(k+1,i+1) * b(i+1)
end do
end do
end function linear
end module dfrank_linear

!****h* /testlinsolve
! PURPOSE
! Test linear solvers.
!****
program testlinsolve
use conf
use f95_lapack
use dfrank_linear

implicit none

! call run_frank_test ()

call run_hilbert_test ()

contains

!****t* testlinsolve/inv_cond
! PURPOSE
! Try to calculate the inverse condition number of a matrix.
!****
subroutine inv_cond (a, rcond, illcond)
real(wp), intent(inout) :: a(:,:)
real(wp), intent(out) :: rcond
logical, optional :: illcond
real(wp) :: s(size(a, 1))

!Do SVD via LAPACK95
call la_gesdd (a, s)

! Condition number of the matrix is the ratio between
! the largest to the smallest singular value. If this
! approaches infinity, the matrix is ill-conditioned,
! i.e. the input vectors were not linearly independent.
! For numerical reasons calculate the inverse of the
! condition number. If this number approaches the machine precision,
! we're screwed.
rcond = minval (s) / maxval (s)
if (rcond < epsilon (rcond)) then
if (present (illcond)) then
illcond = .true.
end if
else
if (present (illcond)) then
illcond = .false.
end if
end if
end subroutine inv_cond

!****t* testlinsolve/hilbert
! PURPOSE
! Make a Hilbert matrix, useful for testing how well linear solvers
! deal with ill-conditioned matrices.
!****
subroutine make_hilbert (a)
real(wp), intent(out) :: a(:,:)
integer :: i, j

forall (i = 1:size (a, 1), j = 1:size (a, 2))
a(i, j) = 1.0_wp / real (i + j - 1, wp)
end forall
end subroutine make_hilbert

!****t* testlinsolve/run_hilbert_test
! PURPOSE
! Test linear solvers using the Hilber matrix.
!****
subroutine run_hilbert_test ()
integer, parameter :: n = 12
real(wp) :: a(n, n), acpy(n,n), rcond, b(n), time1, &
time2, x1(n), b2(n,1), x(n)
logical :: illcond

call make_hilbert (a)
call random_seed ()
call random_number (x) !x is the correct solution
b = matmul (a, x)

b2(:,1) = b
acpy = a

! print *, 'Correct solution is: ', x

call inv_cond (a, rcond, illcond)
if (illcond) then
print *, 'Hilbert matrix of size ', n, ' is ill-conditioned.'
else
print *, 'Hilbert matrix of size ', n, ' is well-conditioned.'
end if
print *, 'rcond: ', rcond
print *, 'Condition number: ', 1.0_wp/rcond

a = acpy
call cpu_time (time1)
call la_gesv (a, b2)
call cpu_time (time2)
print *, 'la_gesv solved system in ', time2 - time1, ' seconds.'
print *, 'Average cumulative error is: ', sum (abs (b2(:,1)-x)) / n
! print *, 'Solution is: ', b2(:,1)

a = acpy
call cpu_time (time1)
x1 = linear (a, b, n-1)
call cpu_time (time2)
print *, 'David Frank solved system in ', time2 - time1, ' seconds.'
print *, 'Average cumulative error is: ', sum (abs (x1-x)) / n
! print *, 'Solution is: ', x1
end subroutine run_hilbert_test

end program testlinsolve

--
Janne Blomqvist

Rich Townsend

2005-06-08, 8:58 pm

Janne Blomqvist wrote:
<snip>

> We see that la_gesv (which uses LU) is more than an order of magnitude
> faster. Well, even LU can't get reasonable result with a Hilbert
> matrix this size, but for size 11 (size 12 happens to be about the
> maximum la_gesv can handle in this case) we see:
>
> % ./testlinsolve
> Hilbert matrix of size 11 is well-conditioned.
> rcond: 1.913805532331239E-15
> Condition number: 522519129089298.4
> la_gesv solved system in 0.000000000000000E+00 seconds.
> Average cumulative error is: 3.605202376467460E-03
> David Frank solved system in 0.000000000000000E+00 seconds.
> Average cumulative error is: 18642906443420.82
>
> While the runtimes are neglible in both cases, Franks code has some
> serious issues with accuracy.

Serious issues? That must be the understatement of the year. We're
talking 17 orders of magnitude here!

David Frank

2005-06-09, 8:57 am

"Janne Blomqvist" <foo@bar.invalid> wrote in message
news:slrndaelos.2g3.foo@vipunen.hut.fi...
> In article <d87dj2$h1$1@scrotar.nss.udel.edu>, Rich Townsend wrote:
>
> Here goes. I split up Davids benchmark, and added another one, as well
> as converting it to f95 so that I could compile it with Lahey.

In reply to Rich,
OK, I did add a timing inside my linear function after the data summation
and find the time taken to invert and
use the invert for coef compute takes less than 1 windows clock tick which
means it takes less than 0.016sec
on my pentium 2.8ghz

Neither of you report what processing my "test" coef with lapack generates
as a
average coef error = sum(abs(coef-test))/(nc+1)
whereas my function's error is 0.000000000000xxx

I can only assume your error is quite a bit larger and you are suppressing
that fact,
if so and since use of invert adds NIL repeat NIL to runtime I fail to see
what there is to complain about.
the use of invert if it is also more accurate for NC < 100 as per my
original claim for the Gauss invert routine
I posted on day 1 of this discussion.

As to stability, I can add a loop to run my linear function with millions of
different random test cases and there wont
be even 1 hiccup where avg coef err isnt contained to a usable result.

Sooo, in the real world, no-one shud expect any problem with their REAL data
sets and #coef < 100
which is what I originally claimed for the Gauss invert algorithm.
1. briefer
2. faster
3. less error

And the winner is: http://home.earthlink.net/~dave_gemini/linear.f90

Greg Lindahl

2005-06-09, 8:57 am

In article <khSpe.1303$eM6.1173@newsread3.news.atl.earthlink.net>,
David Frank <dave_frank@hotmail.com> wrote:

>Sooo, in the real world, no-one shud expect any problem with their REAL data
>sets and #coef < 100
>which is what I originally claimed for the Gauss invert algorithm.
>1. briefer
>2. faster
>3. less error

Just remember, folks, his command of performance issues is as good as
his command of numerical methods. So next time he asks you for a
benchmark number, just say no.

-- greg

Rich Townsend

2005-06-09, 3:58 pm

David Frank wrote:
> "Janne Blomqvist" <foo@bar.invalid> wrote in message
> news:slrndaelos.2g3.foo@vipunen.hut.fi...
>
>
>
> In reply to Rich,
> OK, I did add a timing inside my linear function after the data summation
> and find the time taken to invert and
> use the invert for coef compute takes less than 1 windows clock tick which
> means it takes less than 0.016sec
> on my pentium 2.8ghz
>
> Neither of you report what processing my "test" coef with lapack generates
> as a
> average coef error = sum(abs(coef-test))/(nc+1)
> whereas my function's error is 0.000000000000xxx
>
> I can only assume your error is quite a bit larger and you are suppressing
> that fact,

Spare us the conspiracy theory, David. Here are some results:

linear:
e.g. pentium 2.8ghz, CVF, time = 0.5429
coef recovered = 10.0000 0.5000 0.4900 0.4800 0.4700 0.0100
coef avg error = 0.000000000000125

linear_lapack:
e.g. pentium 2.8ghz, CVF, time = 0.5399
coef recovered = 10.0000 0.5000 0.4900 0.4800 0.4700 0.0100
coef avg error = 0.000000000000059

The actual error obviously changes from execution to execution, but
always seems to involve at most the last 3 digits. Essentially, there's
no difference between the two.

> if so and since use of invert adds NIL repeat NIL to runtime I fail to see
> what there is to complain about.
> the use of invert if it is also more accurate for NC < 100 as per my
> original claim for the Gauss invert routine
> I posted on day 1 of this discussion.

Well, that claim is certainly true. But as has been demonstrated
elsewhere, your matrix inversion approach is slow and unstable. It's
just that your least-squares fitting program is insensitive to these
problems.

>
> As to stability, I can add a loop to run my linear function with millions of
> different random test cases and there wont
> be even 1 hiccup where avg coef err isnt contained to a usable result.

Yes, but your matrix is never nearly deficient.

>
> Sooo, in the real world, no-one shud expect any problem with their REAL data
> sets and #coef < 100
> which is what I originally claimed for the Gauss invert algorithm.
> 1. briefer
> 2. faster
> 3. less error

Let's just review the facts here. My linear_lapack.f90 code takes less
lines of code than your linear.f90 (since the matrix inversion guff is
replaced by a single LAPACK call); as I demonstrated, it is marginally
faster; and there is essentially no difference in the error.

So your points 1-3 are the reverse of the truth.

>
> And the winner is: http://home.earthlink.net/~dave_gemini/linear.f90

Unfortunately not.

cheers,

Rich

Abdul Qat

2005-06-09, 3:58 pm

"David Frank" <dave_frank@hotmail.com> wrote in message
news:khSpe.1303$eM6.1173@newsread3.news.atl.earthlink.net...

[...]

FYI, a Hilbert matrix can be inverted exactly on paper; it's only ill
conditioned in _finite_ precision arithmetic; and it's generally regarded
as a poor test matrix for gauging the performance of linear solvers.

--
You're Welcome,
Gerry T.
______
"Science is the belief in the ignorance of experts." -- Feynman, in The
Pleasure of Finding Things Out.

Sponsored Links