bugGNU Octave - Bugs: bug #60846, sparse matrix elementwise...


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

Submitted by:  John W. Eaton <jwe>
Submitted on:  Mon 28 Jun 2021 08:14:34 PM UTC  
Category:  Octave Function Severity:  4 - Important
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  None Assigned to:  John W. Eaton <jwe>
Originator Name:  jwe Open/Closed:  Open
Release:  dev Operating System:  Any

Add a New Comment (Rich Markup)

You are not logged in

Please log in, so followups can be emailed to you.


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>
Project AdministratorIn charge of this item.


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

Attach Files:

No files currently attached


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by jwe (Submitted the item)

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.




    No changes have been made to this item

    Back to the top

    Powered by Savane 3.6