cholesky decompostion and inverses of upper diagonal matrices 
Author Message
 cholesky decompostion and inverses of upper diagonal matrices

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.



Sun, 25 Nov 2007 01:46:04 GMT  
 cholesky decompostion and inverses of upper diagonal matrices
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



Sun, 25 Nov 2007 02:32:29 GMT  
 cholesky decompostion and inverses of upper diagonal matrices

Quote:

> 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



Sun, 25 Nov 2007 18:38:37 GMT  
 cholesky decompostion and inverses of upper diagonal matrices


<snip>

Quote:

> 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



Sun, 25 Nov 2007 21:28:36 GMT  
 cholesky decompostion and inverses of upper diagonal matrices

Quote:

> 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



Mon, 26 Nov 2007 01:08:53 GMT  
 cholesky decompostion and inverses of upper diagonal matrices

Quote:



> <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

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

:)

cheers,

Rich



Mon, 26 Nov 2007 01:45:40 GMT  
 cholesky decompostion and inverses of upper diagonal matrices


Quote:



>> <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

> 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



Mon, 26 Nov 2007 02:01:34 GMT  
 cholesky decompostion and inverses of upper diagonal matrices

Quote:






>>><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

>>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

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



Mon, 26 Nov 2007 02:32:17 GMT  
 cholesky decompostion and inverses of upper diagonal matrices

Quote:

> 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



Mon, 26 Nov 2007 04:37:48 GMT  
 cholesky decompostion and inverses of upper diagonal matrices
<snip>

Quote:
> 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!


Mon, 26 Nov 2007 05:04:11 GMT  
 cholesky decompostion and inverses of upper diagonal matrices


Quote:

>> 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.

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



Mon, 26 Nov 2007 15:33:36 GMT  
 cholesky decompostion and inverses of upper diagonal matrices

Quote:

>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



Mon, 26 Nov 2007 15:50:44 GMT  
 cholesky decompostion and inverses of upper diagonal matrices

Quote:




>>>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.

> 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 {*filter*} 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.

Quote:
> 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.

Quote:

> 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.

Quote:

> 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.

Quote:

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

Unfortunately not.

cheers,

Rich



Mon, 26 Nov 2007 21:00:54 GMT  
 cholesky decompostion and inverses of upper diagonal matrices


[...]

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.



Mon, 26 Nov 2007 23:30:23 GMT  
 
 [ 14 post ] 

 Relevant Pages 

1. Cholesky inverse of a covariance matrix

2. Question: Tri-diagonal matrices in APL

3. LOGO-L> Diagonal Matrix

4. Challenge 4 (matrix diagonal)

5. Eigenvectors of tri-diagonal matrix

6. Eigenvectors for tri-diagonal matrix

7. Initialize a Diagonal Matrix

8. diagonal matrix partial sums

9. APL matrix inverse

10. Matrix Inverse in APL2

11. Matrix#inverse broken?

12. Finding the inverse of a matrix.

 

 
Powered by phpBB® Forum Software