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

 
 

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

Submitter:  John W. Eaton <jwe>
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
   

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

 

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

Attach Files:
   
   
Comment:
   

No files currently attached

 

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.

     

    Follow 2 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2024-01-28 rik5 StatusNone Confirmed
    2024-01-28 rik5 Dependencies- Depends on bugs #57909

    Back to the top

    Powered by Savane 3.13-b921.
    Corresponding source code