[Octave-bug-tracker] [bug #60846] sparse matrix elementwise exponentiation can produce incorrect results

[John W. Eaton]

URL:
  <https://savannah.gnu.org/bugs/?60846>

                 Summary: sparse matrix elementwise exponentiation can produce
incorrect results
                 Project: GNU Octave
            Submitted by: 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: jwe
         Originator Name: jwe
        Originator Email: 
             Open/Closed: Open
                 Release: dev
         Discussion Lock: Any
        Operating System: Any

    _______________________________________________________

Details:

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.







    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?60846>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/