# bug #60846: sparse matrix elementwise exponentiation can produce incorrect results

 Submitter: John W. Eaton Submitted: Mon 28 Jun 2021 08:14:34 PM UTC Category: Octave Function Severity: 4 - Important Priority: 5 - Normal Item Group: Incorrect Result Status: Confirmed Assigned to: jwe 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 12:45:57 AM UTC, comment #1:

Added bug #57909 as a dependency.  There is already a patch submitted there.  Either this approach, or the one at the other bug report should be used.

Rik <rik5>
Mon 28 Jun 2021 08:14:34 PM UTC, original submission:

I see an incorrect result for the following expression:

sparse(0) .^ sparse(1)  ==>  sparse(1) but should should be sparse(0)

While investigating this error, I also found the following incorrect result:

octave> [0,1;0,2] .^ [0,0;1,2]  %% full is correct
ans =

1   1
0   4

%% Using full in the expression below so comparison is easier:
octave> full (sparse([0,1;0,2]) .^ sparse([0,0;1,2]))
ans =

1   1
1   4

Since the 1x1 sparse matrix objects are not converted to scalars, the following function in sparse-xpow.cc is used:

elem_xpow (const SparseMatrix& a, const SparseMatrix& b)

and we ultimately end up in the following block at the end of the function:

SparseMatrix result (nr, nc, 1.0);

for (octave_idx_type j = 0; j < nc; j++)
{
for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
{
octave_quit ();
result.xelem (a.ridx (i), j) = std::pow (a.data (i),
b(a.ridx (i), j));
}
}
result.maybe_compress (true);

If I understand correctly, this loop is only operating on the non-zero elements of the LHS of the expression and not accounting for the possibility of a 0 element on the LHS being raised to a non-zero power.  It seems we need something more like the loops used in Sparse-op-defs.h that use the following pattern:

r = SparseBoolMatrix (m1_nr, m1_nc, true);              \
for (octave_idx_type j = 0; j < m1_nc; j++)             \
{                                                     \
octave_idx_type i1 = m1.cidx (j);                   \
octave_idx_type e1 = m1.cidx (j+1);                 \
octave_idx_type i2 = m2.cidx (j);                   \
octave_idx_type e2 = m2.cidx (j+1);                 \
while (i1 < e1 || i2 < e2)                          \
{                                                 \
if (i1 == e1 || (i2 < e2 && m1.ridx (i1) > m2.ridx (i2))) \
{                                             \
if (! (Z1 OP m2.data (i2)))                 \
r.data (m2.ridx (i2) + j * m1_nr) = false; \
i2++;                                       \
}                                             \
else if (i2 == e2 || m1.ridx (i1) < m2.ridx (i2)) \
{                                             \
if (! (m1.data (i1) OP Z2))                 \
r.data (m1.ridx (i1) + j * m1_nr) = false; \
i1++;                                       \
}                                             \
else                                            \
{                                             \
if (! (m1.data (i1) OP m2.data (i2)))       \
r.data (m1.ridx (i1) + j * m1_nr) = false; \
i1++;                                       \
i2++;                                       \
}                                             \
}                                                 \
}                                                     \
r.maybe_compress (true);

This way, we will touch all the cases where either the LHS or the RHS (or both) have non-zero elements.

Also in the SPARSE_SMSM_CMP_OP macro in Sparse-op-defs.h, I see that there are two branches, one for the case when the comparison operator returns TRUE for LHS and RHS both zero (that's shown above) and another for when it returns FALSE.  It's not clear to me that having two branches is really necessary.

And, like the recent discussion for mpower in bug #60786, it would probably be best to convert the code in Sparse-op-defs.h to use templates.  Likewise for the functions in sparse-xpow.cc if possible, though I'm not sure we have as much duplication there as in xpow.cc since we don't have the duplication for single-precision sparse matrices.

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:

No files currently attached

## Dependencies

Depends on the following items

Digest:
bug dependencies.

Items that depend on this one: None found

Carbon-Copy List
• -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.