# bug #57909: incorrect result for sparse matrix exponentiation

 Submitter: John W. Eaton Submitted: Thu 27 Feb 2020 05:56:23 PM UTC Category: Octave Function Severity: 4 - Important Priority: 7 - High Item Group: Incorrect Result Status: Patch Submitted Assigned to: None Originator Name: jwe Open/Closed:  * Open Release:  * dev Operating System:  * Any Fixed Release: None Planned Release: None
* Mandatory Fields

Add a New Comment Rich Markup

## Discussion

Sun 28 Jan 2024 01:25:30 AM UTC, comment #9:

This is correct

c = as .^ 1
c = Compressed Column Sparse (rows = 1, cols = 1, nnz = 0 [0%])

Check for full matrices

octave:3> 0 .^ 1
ans = 0

And sparse display only reports non-zero entries.  For example,

octave:4> sparse ([0])
ans = Compressed Column Sparse (rows = 1, cols = 1, nnz = 0 [0%])

Size is 1x1, and number of non-zero entries is zero.

Rik <rik5>
Sun 28 Jan 2024 01:03:21 AM UTC, comment #8:

But that looks like a variable display issue (a different bug?)

octave:17> a= zeros(2)
a =

0   0
0   0

octave:18> as = sparse (a)
as =

Compressed Column Sparse (rows = 2, cols = 2, nnz = 0 [0%])

octave:19> aa = full(as)
aa =

0   0
0   0

octave:20> as(1)
ans = Compressed Column Sparse (rows = 1, cols = 1, nnz = 0 [0%])

octave:21> as(1)+1
ans = 1
octave:22> 1+ as(1)
ans = 1
octave:23> whos as
Variables visible from the current scope:

variables in scope: top scope

Attr   Name        Size                     Bytes  Class
====   ====        ====                     =====  =====
s     as          2x2                         40  double

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Sun 28 Jan 2024 12:57:19 AM UTC, comment #7:

octave:9> c = as .^ 1
c = Compressed Column Sparse (rows = 1, cols = 1, nnz = 0 [0%])

Dmitri
--

Dmitri A. Sergatskov <dasergatskov>
Sun 28 Jan 2024 12:43:38 AM UTC, comment #6:

Not sure why this wasn't reviewed four years ago, but this is still a real problem.  Minimum example below.

octave:1> a = 0
a = 0
octave:2> b = 2
b = 2
octave:3> a .^ b
ans = 0
octave:4> as = sparse (a)
as = Compressed Column Sparse (rows = 1, cols = 1, nnz = 0 [0%])

octave:5> bs = sparse (b)
bs = Compressed Column Sparse (rows = 1, cols = 1, nnz = 1 [100%])

(1, 1) -> 2
octave:6> as .^ bs
ans = Compressed Column Sparse (rows = 1, cols = 1, nnz = 1 [100%])

(1, 1) -> 1

Rik <rik5>
Mon 12 Sep 2022 01:43:30 AM UTC, comment #5:

arungiridhar is working on sparse-xpow to improve speed

https://octave.discourse.group/t/why-is-sparse-matrix-exponentiation-slower-than-iterated-multiplication/3216

Should jwe's patch and bug #60846 be included their investigation?

Anonymous
Fri 28 Feb 2020 03:39:30 PM UTC, comment #4:

I'm attaching a patch to fix the functions in sparse-xpow.cc.  I'll need to also add tests.  I'd be glad to see any improvements in either the algorithm or the use of templates for this function.

(file #48517)

John W. Eaton <jwe>
Thu 27 Feb 2020 07:30:45 PM UTC, comment #3:

Matlab 2019a:

>> a = sparse (4, 6);
b = sparse (4, 6);
a(2,2) = 1;
b(3,3) = 2;
a(3,4) = 3; b(2,4) = 4;
a(2,5) = 5; b(3,5) = 6;
a(4,6) = 7; b(4,6) = 8;
>> full (a .^ b)
ans =
1           1           1           1           1           1
1           1           1           0           1           1
1           1           0           1           0           1
1           1           1           1           1     5764801
>> full (a) .^ full (b)
ans =
1           1           1           1           1           1
1           1           1           0           1           1
1           1           0           1           0           1
1           1           1           1           1     5764801

>> a = sparse (0)
a =
All zero sparse: 1×1
>> b = sparse (2)
b =
(1,1)        2
>> a.^b
ans =
All zero sparse: 1×1
>> 0.^2
ans =
0

Nicholas Jankowski <nrjank>
Thu 27 Feb 2020 07:26:28 PM UTC, comment #2:

I think the attached function does the right indexing when there is no need to convert to complex.  A proper replacement for the elem_xpow function in sparse-xpow.cc also needs to handle the complex case.  And since the indexing should be the same for the real and complex cases, maybe it should be done with a template rather than duplicating the code.

(file #48511)

John W. Eaton <jwe>
Thu 27 Feb 2020 06:46:58 PM UTC, comment #1:

Confirmed.  This doesn't appear correct.  For a more vivid example, just take one of the failing entries such as a(3,3) => 0 and b(3,3) = 2.

octave:14> a = sparse (0)
a = Compressed Column Sparse (rows = 1, cols = 1, nnz = 0 [0%])

octave:15> b = sparse (2)
b = Compressed Column Sparse (rows = 1, cols = 1, nnz = 1 [100%])

(1, 1) -> 2
octave:17> a.^b
ans = Compressed Column Sparse (rows = 1, cols = 1, nnz = 1 [100%])

(1, 1) -> 1

But, when calculated in a regular manner

octave:20> 0 .^ 2
ans = 0

Rik <rik5>
Thu 27 Feb 2020 05:56:23 PM UTC, original submission:

The element-by-element sparse exponentiation functions miss the case when the LHS is 0 and the RHS is non-zero and the result should be zero.  In the following, I believe the results should agree for both sparse and full, and that the sparse result is wrong:

a = sparse (4, 6);
b = sparse (4, 6);
a(2,2) = 1;
b(3,3) = 2;
a(3,4) = 3; b(2,4) = 4;
a(2,5) = 5; b(3,5) = 6;
a(4,6) = 7; b(4,6) = 8;
full (a .^ b) =>

1         1         1         1         1         1
1         1         1         1         1         1
1         1         1         1         1         1
1         1         1         1         1   5764801

full (a) .^ full (b) =>

1         1         1         1         1         1
1         1         1         0         1         1
1         1         0         1         0         1
1         1         1         1         1   5764801

John W. Eaton <jwe>

## Attached Files

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

Comment:

Attached Files
file #48517:  sparse-pow-diffs.txt added by jwe (7KiB - text/plain)
file #48511:  sparsexpow.cc added by jwe (3KiB - text/x-c++src)

## Dependencies

Depends on the following items: None found

Items that depend on this one

Digest:
bug dependencies.

Carbon-Copy List
• -email is unavailable- added by dasergatskov (Posted a comment)
• -email is unavailable- added by nrjank (Posted a comment)
• -email is unavailable- added by rik5 (Posted a comment)
• -email is unavailable- added by jwe (Submitted the item)
•

There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

Only group members can vote.