sparse symmetric matrices

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

sparse symmetric matrices

Bill Greene
I'm trying to solve a linear system with a sparse,
symmetric matrix and have some questions I'm hoping
someone can help me with. I'm using Octave-3.2.4,
by the way.

I'm calling chol() directly and providing only entries
in the upper triangle of the matrix. This works fine
even though the Octave doc does not say specifically that
only the upper triangle is needed.

However, if I first call amd() to order the matrix, I get
a permutation that puts the entries in both upper and lower
triangles so chol() fails.

My questions are:

1) Is there a fundamentally better way to deal with this type
    of system (I strongly prefer storing only 1/2 the matrix)?
2) Is there a method (with good performance) that will
    rearrange my permuted sparse matrix to move all the lower
    entries to the upper triangle?

Here is a code sample that works:
K = [[0,0,2]; [0,1,-2]; [0,4, -1];
       [1,1,3]; [1,2,-2];
       [2,2,5]; [2,3,-3]; [3,3,10]; [3,4,4]; [4,4,10]];
  for i = 1:2
    K(:,i) += 1;
  end
  ks = sparse(K(:,1), K(:,2), K(:,3));
  b = [0,1,0,0,0];
  lt = chol(ks);
  u = lt \ (lt'\b')

But if I replace the relevant lines with:
p = amd(ks)
lt = chol(ks(p,p));
x = lt \ (lt'\b(p)');
u = x(p);

I run into the problem described above.

Thanks.

Bill Greene

_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: sparse symmetric matrices

Bill Greene-2
I guess I should have mentioned that I did try:

u = ks \ b'

This returns an incorrect solution, u.
I also tried to let chol() compute a permutation
[lt, err, p] = chol(ks, 'vector');
and this results in an incorrect solution, also.

I'm guessing all of these errors are related: I'm
creating only the upper triangle and, when permuted,
entries end up in both triangles.

I'm hoping someone has a good strategy for dealing with
this.

Bill


On 6/11/2011 1:48 PM, David Bateman wrote:

> On 06/11/2011 07:11 PM, Bill Greene wrote:
>> I'm trying to solve a linear system with a sparse,
>> symmetric matrix and have some questions I'm hoping
>> someone can help me with. I'm using Octave-3.2.4,
>> by the way.
>>
>> I'm calling chol() directly and providing only entries
>> in the upper triangle of the matrix. This works fine
>> even though the Octave doc does not say specifically that
>> only the upper triangle is needed.
>>
>> However, if I first call amd() to order the matrix, I get
>> a permutation that puts the entries in both upper and lower
>> triangles so chol() fails.
>>
>> My questions are:
>>
>> 1) Is there a fundamentally better way to deal with this type
>>     of system (I strongly prefer storing only 1/2 the matrix)?
>> 2) Is there a method (with good performance) that will
>>     rearrange my permuted sparse matrix to move all the lower
>>     entries to the upper triangle?
>>
>> Here is a code sample that works:
>> K = [[0,0,2]; [0,1,-2]; [0,4, -1];
>>        [1,1,3]; [1,2,-2];
>>        [2,2,5]; [2,3,-3]; [3,3,10]; [3,4,4]; [4,4,10]];
>>   for i = 1:2
>>     K(:,i) += 1;
>>   end
>>   ks = sparse(K(:,1), K(:,2), K(:,3));
>>   b = [0,1,0,0,0];
>>   lt = chol(ks);
>>   u = lt \ (lt'\b')
>>
>> But if I replace the relevant lines with:
>> p = amd(ks)
>> lt = chol(ks(p,p));
>> x = lt \ (lt'\b(p)');
>> u = x(p);
>>
>> I run into the problem described above.
>>
>> Thanks.
>>
>> Bill Greene
> Why not just call
>
> x = ks \ b;
>
> and let Octave handle the recognition that the matrix is positive
> definite, the permutation, the cholesky and the back substitution. If
> the time needed to recognize that the matrix is positive definite is too
> long. you can call matrix_type on ks before solving the problem to
> bypass the automatic recognition of the matrix type
>
> D.
>
>
>

_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: sparse symmetric matrices

David Bateman
On 06/11/2011 09:38 PM, Bill Greene wrote:

> I guess I should have mentioned that I did try:
>
> u = ks \ b'
>
> This returns an incorrect solution, u.
> I also tried to let chol() compute a permutation
> [lt, err, p] = chol(ks, 'vector');
> and this results in an incorrect solution, also.
>
> I'm guessing all of these errors are related: I'm
> creating only the upper triangle and, when permuted,
> entries end up in both triangles.
>
> I'm hoping someone has a good strategy for dealing with
> this.
>
> Bill
>

Because Octave is recognising that the matrix is upper triangular and
treating it as such, even though the choloesky only uses the upper
triangular part of the matrix. I only see one solution

Tell Octave that the matrix is positive definite rather than upper
triangular. cf

octave:1> n = 1000; d = 0.1; a = sprand(n,n,d); s0 = matrix_type(a *
a',"positive definite"); x = rand(n,1); y0 = s0 \ x; s1 =
matrix_type(triu(s0), "positive definite"); y1 = s1 \ x; norm(y1 - y0)
ans = 0


D.
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave