bugGNU Octave - Bugs: bug #60784, Inconsistent behavior for boolean...

 
 

bug #60784: Inconsistent behavior for boolean matrix types with matrix functions

Submitter:  None
Submitted:  Wed 16 Jun 2021 01:36:01 AM UTC
   
 
Category:  Octave Function Severity:  1 - Wish
Priority:  1 - Later Item Group:  Unexpected Error or Warning
Status:  Need Info Assigned to:  None
Originator Name:  Originator Email:  -email is unavailable-
Open/Closed:  * Open Release:  * dev
Operating System:  * Any Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Tue 22 Jun 2021 05:38:36 PM UTC, comment #14: 

Automatic conversion of integer types to double is probably not desirable because of the saturation semantics of other integer operations:


uint8 (200) + uint8 (200) ==> uint8 (255), not double (400)


John W. Eaton <jwe>
Group administrator
Tue 22 Jun 2021 04:37:04 PM UTC, comment #13: 

To summarize the discussion so far:

1. The following functions automatically convert logical and integer types to double in Octave but not Matlab: det, lu, svd, inv.

2. The eig function consistently errors for anything other than single and double, in both Octave and Matlab.

Possible courses of action to make the functions eig, det, lu, svd and inv all consistent within Octave:

1. If we want to emulate Matlab, throw an error for logical or integer types with det, lu, svd or inv. This will make those functions behave like eig.

2. If we want to be a superset of Matlab, make eig behave like det, inv, lu or svd in accepting logical or integer types. If Octave does silent conversion to double in cases that Matlab doesn't, it could emit a warning so that the user can either suppress that warning ID or do an explicit cast or something different.

Once we agree on one of the above, we have another choice: do we make the integer types automatically convert to double like boolean, or not?

Anonymous
Thu 17 Jun 2021 02:04:12 AM UTC, comment #12: 

Thanks for the test (comment #11) as more solid basis for discussion.  In Matlab R2021a


Type logical:
    Silent: add, subtract, matrix multiply, matrix ^ scalar, scalar ^ matrix
    Error: inv, det, lu, eig, svd

Types int8, int16, int32, int64, uint8, uint16, uint32, uint64:
    Silent: add, subtract
    Error: matrix multiply, matrix ^ scalar, scalar ^ matrix, inv, det, lu, eig, svd

## Hint "matrix multiply A*B"
## Error using  *
## MTIMES (*) is not fully supported for integer classes. At least one argument must be scalar.

Types single, double:
    All functions work, no errors or warnings.



Kai Torben Ohlhus <siko1056>
Group Member
Wed 16 Jun 2021 10:36:29 PM UTC, comment #11: 

From curiosity, I did a cross product of multiple types with multiple functions using the code at the bottom of this comment. Here are the results:


Type logical:
    Silent: add, subtract, matrix multiply, matrix ^ scalar, scalar ^ matrix, det, lu, svd
    Warning: inv (singular to machine precision)
    Error: eig

Types int8, int16, int32, int64, uint8, uint16, uint32, uint64:
    Silent: add, subtract, inv, det, lu, svd
    Warning:
    Error: matrix multiply, matrix ^ scalar, scalar ^ matrix, eig

Types single, double:
    All functions work, no errors or warnings.


Code to test:

origA = rand(50); # type double
origB = rand(50); # type double

types = ["logical"; "int8"; "int16"; "int32"; "int64"; "uint8"; "uint16"; "uint32"; "uint64"; "single"; "double"];

funcs = ["A+B"; "A-B"; "A*B"; "A^1.6"; "(-1.6)^A"; "inv(A)"; "det(A)"; "lu(A)"; "eig(A)"; "svd(A)"];

for i = 1:size(types,1)
  for j = 1:size(funcs,1)
    str = sprintf("A = %s(origA); B = %s(origB); C = %s;",...
                  types(i,:), types(i,:), funcs(j,:));
    disp(str)
    try
      eval(str)
    catch
      disp("error for the last operation")
    end
    disp(' ')
  end
end


Anonymous
Wed 16 Jun 2021 09:57:39 PM UTC, comment #10: 

While this discussion is going on, here is some more material to discuss:


A = (rand(1000) > 0.5); # A is boolean
B = A^2;                # boolean^2 converts to double
B = single(A)^2;        # B is type single
B = int8(A)^2;          # error
B = int16(A)^2;         # error
B = int32(A)^2;         # error


Now why is it OK to silently convert a boolean type to double but not the integer types? Opening up that possibility could create problems with memory usage though, so it's not clear-cut.

Based on comment #1 and comment #9, it might be better to require that all matrix functions require single or double, and any type conversions need to be done by the user explicitly. The price would be user inconvenience though, so again the tradeoff is not clear.

Anonymous
Wed 16 Jun 2021 03:49:17 PM UTC, comment #9: 

Prior to this change (from 2008!)

http://hg.savannah.gnu.org/hgweb/octave/diff/18c4ded8612a/src/DLD-FUNCTIONS/eig.cc

it looks like eig would have worked on integer and logical matrices.  But then the test was added to error early if the argument was not single or double, then later that test was changed to use isfloat.

In any case, can someone test what Matlab actually does for the following?


eig (rand (3) > 0.5)
eig (int8 (rand (3) * 255))
eig (int16 (rand (3) * 65535))


I'm not saying we can't support this feature but I'm wondering whether Matlab actually restricts the inputs to double and single as the docs say, or if it silently converts other numeric types to double?

If it errors for types other than single and double and we accept  and convert other numeric types, then how long will it be until someone submits a bug report about a missed (instead of unexpected) warning or error?  LOL.

John W. Eaton <jwe>
Group administrator
Wed 16 Jun 2021 04:45:01 AM UTC, comment #8: 

"Basic mem." is the main memory occupied by the operating system and after starting Octave.
"def. A" is the main memory usage after defining the Matrix A.
"det (A)" is the main memory usage while computing det(A).
"det (double (A))" like above, only done for the "logical" case.

You can use a tool like htop and observe on a second terminal screen (not as convenient as tic-toc).

The problem is, there are only "single" and "double" implementations available.  Everything not "single" from the start gets eventually casted to "double" under the hood.

This is what I mean by "no surprises".  The user has to decide which of the two available LAPACK implementations "single" or "double" is suited for his problem.

Kai Torben Ohlhus <siko1056>
Group Member
Wed 16 Jun 2021 04:39:42 AM UTC, comment #7: 

I am not sure how do you measure memory. Also what is "def. A"?

My point was that "optimal" casting depends on memory vs speed considerations (single is faster than int16, but takes more memory).

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Wed 16 Jun 2021 04:26:13 AM UTC, comment #6: 

Adding a timing might have obfuscated the point I tried to make, sorry for the confusion 😓

On a machine with 1 TB of main memory:


A = rand (2e4);
B = logical (A);
C = single (B);
D = int16(B);

>> whos
Variables visible from the current scope:

variables in scope: top scope

   Attr Name        Size                     Bytes  Class
   ==== ====        ====                     =====  =====
        A       20000x20000             3200000000  double
        B       20000x20000              400000000  logical
        C       20000x20000             1600000000  single
        D       20000x20000              800000000  int16
        ans         1x1                          8  double

Total is 1600000001 elements using 6000000008 bytes


tic; det (A); toc           # Elapsed time is 24.2933 seconds.
tic; det (B); toc           # Elapsed time is 14.8445 seconds.
tic; det (double (B)); toc  # Elapsed time is 15.7939 seconds.
tic; det (C); toc           # Elapsed time is 8.57522 seconds.
tic; det (D); toc           # Elapsed time is 14.6047 seconds.


Thus roughly confirmed your timing observations.  The difference with det(A) to det(B) I explain with the "nicer" structure of B (only ones) the algorithm can exploit (DGETRF, LAPACK LU-factorization, summing up the diagonal)
https://hg.savannah.gnu.org/hgweb/octave/file/84d7e1ace446/liboctave/array/dMatrix.cc#l862

My point is memory consumption during the algorithm execution.  The timings themselves are of minor interest and I introduced them only to "ensure" that the same "branches" of the algorithm are taken (not a smart decision apparently).

Now when monitoring only the memory


##                 Basic mem.   def. A   det (A)                  det (double (A))
A = rand (2e4);  #   2.1 GB     5.1 GB   8.0 GB (10.0 GB peak)
A = logical (A); #   2.1 GB     2.5 GB   8.4 GB (11.4 GB peak)    8.4 GB (11.4 GB peak)
A = single (A);  #   2.1 GB     3.6 GB   5.1 GB ( 6.1 GB peak)
A = int16 (A);   #   2.1 GB     2.8 GB   8.8 GB (11.8 GB peak)


Can you make similar observations @dasergatskov?

Regarding comment #5: I think you prefer user convenience then?

Kai Torben Ohlhus <siko1056>
Group Member
Wed 16 Jun 2021 04:16:22 AM UTC, comment #5: 

RE comment #3:

If you look at det.cc or other functions like this, I think the general pattern is to test for specific things like diagonal or permutation matrices or single types if they can be handled as special cases (for example, single precision lapack functions) otherwise, we end up at tests for real or complex and conversion to double complex or double for everything else.

John W. Eaton <jwe>
Group administrator
Wed 16 Jun 2021 04:12:47 AM UTC, comment #4: 

Looking at libinterp/corefcn/eig.cc, it contains the test


  if (! arg_a.isfloat ())
    err_wrong_type_arg ("eig", arg_a);


and isfloat is only true for single or double objects.  I'll bet that's a very old check?  Anyway, if it is removed, eig will work because the rest of the function allows conversions if possible.

John W. Eaton <jwe>
Group administrator
Wed 16 Jun 2021 03:12:48 AM UTC, comment #3: 

And perhaps a somewhat surprising thing is:


octave:11> D=int16(B);
octave:12> whos
Variables visible from the current scope:

variables in scope: top scope

   Attr Name        Size                     Bytes  Class
   ==== ====        ====                     =====  =====
        A       20000x20000             3200000000  double
        B       20000x20000              400000000  logical
        C       20000x20000             1600000000  single
        D       20000x20000              800000000  int16
        ans         1x1                          8  double

Total is 1600000001 elements using 6000000008 bytes

octave:13> tic; det (D), toc
ans = 0
Elapsed time is 14.9294 seconds.
octave:14>


It appears that int16 gets promoted to double, while single apparently not.


Dmitri.
--
 



Dmitri A. Sergatskov <dasergatskov>
Wed 16 Jun 2021 03:00:51 AM UTC, comment #2: 

You must be hitting swap or something:


octave:1> A = rand (2e4);
octave:2> B = logical (A);
octave:3> C = single(B);
octave:4> whos
Variables visible from the current scope:

variables in scope: top scope

   Attr Name        Size                     Bytes  Class
   ==== ====        ====                     =====  =====
        A       20000x20000             3200000000  double
        B       20000x20000              400000000  logical
        C       20000x20000             1600000000  single

Total is 1200000000 elements using 5200000000 bytes

octave:5>  tic; det (A); toc
Elapsed time is 16.9533 seconds.
octave:6>  tic; det (B), toc
ans = 0
Elapsed time is 14.8917 seconds.
octave:7> tic; det (double (B)), toc
ans = 0
Elapsed time is 14.948 seconds.
octave:8> tic; det (C), toc
ans = 0
Elapsed time is 7.50517 seconds.
octave:9>  tic; det (A), toc
ans = Inf
Elapsed time is 16.9971 seconds.
octave:10>


Dmitri A. Sergatskov <dasergatskov>
Wed 16 Jun 2021 02:39:53 AM UTC, comment #1: 

Thanks for the report.  I agree, that inconsistency can be annoying.  In this case I tend to support the current Matlab R2021a implementation


>> A = logical (rand (4));
>> det (A)
Check for incorrect argument data type or missing argument in call to function 'det'.

>> eig (A)
Error using eig
Invalid data type. Input matrix must be double or single.


(same for inv(A), svd(A), and lu(A)) to strictly forbid calling those function with non-single or non-double input.

It is more for user convenience, but making a type cast explicit helps in understanding memory issues and avoids ugly surprises.  For example


>> A = rand (2e4);   ##   3 GB
>> B = logical (A);  ## 380 MB

>> tic; det (A), toc
ans = -Inf
Elapsed time is 223.995 seconds.
>> tic; det (B), toc
ans = 0
Elapsed time is 136.939 seconds.
tic; det (double (B)), toc
ans = 0
Elapsed time is 135.042 seconds.


An unaware programmer might think if "det (B)" is permitted, there must be a memory efficient implementation for "logical" matrices in Octave.

However, running the example and observing a memory monitor shows a slightly bigger memory blowup (as the original 380 MB matrix B has to be stored on top of the internal conversion to a full double matrix).  Thus one gets an even worse memory usage, even though I did everything to save memory using a logical array 😓

In my opinion, the discussion boils down to convenience vs. no surprises.  As said before, I prefer no surprises.

How do other maintainer think about it? 🙂  Is it worth ensuring consistency and if yes, convenience or no surprises (Matlab compatibility)?

Kai Torben Ohlhus <siko1056>
Group Member
Wed 16 Jun 2021 01:36:01 AM UTC, original submission:  

If A is a boolean matrix, then most matrix operations like addition, multiplication, exponentiation, det(A), inv(A), svd(A), lu(A) etc are OK (the matrix gets converted to numerical type first; certain functions only apply to square matrices of course) but other operations like eig(A) cause type errors that they won't work on boolean matrices.

Is there a rationale for which functions require a strictly floating point numeric type and which ones will do the conversion automatically?

This behavior was encountered when calculating the eigenvalues of a graph's adjacency matrix which was a 0/1 logical type. It was worked around with "eig(double(A))" but it would be nice if all matrix functions behaved the same in whether they accept boolean input or not.

Anonymous

 

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

Attach Files:
   
   
Comment:
   

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 (Posted a comment)
  • -email is unavailable- added by dasergatskov (Posted a comment)
  • -email is unavailable- added by siko1056 (Posted a comment)
  •  

    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 3 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2021-06-16 siko1056 Severity3 - Normal 1 - Wish
        Priority5 - Normal 1 - Later
        StatusNone Need Info

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code