# bug #47415: out of memory negating a permutation matrix

 Submitter: Carlo de Falco Submitted: Tue 15 Mar 2016 09:23:46 AM UTC Category: Octave Function Severity: 3 - Normal Priority: 5 - Normal Item Group: Regression Status: None Assigned to: None Originator Name: Open/Closed:  * Open Release:  * dev Operating System:  * Any Fixed Release: None Planned Release: None
* Mandatory Fields

## Discussion

Fri 16 Sep 2016 11:34:55 PM UTC, comment #32:

It appears that one possibility for maintaining Matlab compatibility would be to have the creation of Permutation matrices depend on the full/sparse quality of the input matrix.

For lu, the motivating example for this bug report, if the input is Sparse we output a Permutation matrix.  This permutation matrix can either convert to sparse when any unsupported operation is asked of it.  Or we could at least implement the unary minus operator (-) for Permutation matrices since this would be quite easy.  This matches Matlab which outputs a Sparse matrix which has no concept of -0.

If lu had a full matrix input, then presumably the problem was not so large to begin with that the outputs L, U, and P can not also be represented by full matrices.  In this case, we just waste the memory and represent P as a full matrix the same as Matlab.

Does that sound right?

Rik <rik5> Thu 15 Sep 2016 03:31:28 PM UTC, comment #31:

issparse(P) returns zero.

in the workspace, P is listed as a 100x100 double.

Nicholas Jankowski <nrjank> Wed 14 Sep 2016 04:42:07 PM UTC, comment #30:

@Nicholas: Could you try the following code sample under Matlab?

A = sprandn (100, 100, .02);
Af = full (A);
[L, U, P] = lu (Af);
issparse (P)

Rik <rik5> Mon 29 Aug 2016 07:02:52 PM UTC, comment #29:

thought I'd add some info but both :

>> [L, U, P, Q, R] = lu (Af);

and

>> [L, U, P, Q] = lu (Af);

produce the error

Error using lu
Too many output arguments.

Nicholas Jankowski <nrjank> Sun 27 Mar 2016 10:06:49 PM UTC, comment #28:

I think Carlo is on to something here.  At least when the input matrix is sparse, the output matrix P should also be sparse (or behave like a sparse matrix) for compatibility.

According to Octave's documentation,

If A is full subroutines from LAPACK are used and if A is sparse then UMFPACK is used.

So there may still be a compatibility issue if Matlab returns a full permutation matrix P when the input is a full matrix.  What does Matlab do for this

A = sprandn (100, 100, .02);
Af = full (A);
[L, U, P, Q, R] = lu (Af);
issparse (P)

Also, if I understand things correctly it shouldn't return the Q and R values as those are only for sparse matrices.

Rik <rik5> Sun 27 Mar 2016 09:32:12 PM UTC, comment #27:

jwe,

also from the point of view of maintaining compatibility
w.r.t. to -0 and -Inf, the conversion to sparse appears to be more correct in this case (factorization reorderings) than a conversion to full.

In Matlab 2015b I see the following:

>> A = sprandn (100, 100, .02);
>> [L, U, P, Q, R] = lu (A);
>> 1./(-P(1:2,1:2))

ans =

Inf   Inf
Inf   Inf

>>

In Octave 4.1 the same code returns

>> A = sprandn (100, 100, .02);
>> [L, U, P, Q, R] = lu (A);
>> 1./(-P(1:2,1:2))
ans =

-Inf  -Inf
-Inf  -Inf

while the conversion to sparse would lead to the following

>> 1./(-sparse(P(1:2,1:2)))
ans =

Inf   Inf
Inf   Inf

honestly the idea of using n^2 instead of ~n memory to store
a lot of '-0's sounds a bit exotic to me, if we really cared
so much about -0 we could add a signbit property to SparseMatrix
as well, but is that really worth the effort??

Carlo de Falco <cdf> Sun 27 Mar 2016 08:56:58 PM UTC, comment #26:

Lachlan,

You are right that diagonal matrices are not
a specialization of sparse, sorry.

As for removing permutation matrices I really hope we
don't choose to go that way!

The reason for adding that specialization, I guess,
was exactly the use case that brought me to discover
this bug, i.e., storing the factorization (with reordering)
of a sparse matrix for later use:

[L, U, P, Q, R] = lu (A);
for ii = 1 : large_number
res = compute_res (ii);
x = Q * (U \ (L \ (P * (R \ res))));
end

which is quite common when solving nonlinear or time dependent PDEs.

Actually I think this (factorization algorithms returning a reordering)
is almost the only example where the average user will get
a permutation matrix.

Running the same factorization in Matlab I get the following

>> A = sprandn (100, 100, .02);
>> [L, U, P, Q, R] = lu (A);
>> issparse (P)
ans =
1
>> issparse (Q)
ans =
1

So Matlab produces directly sparse matrices for the reorderings,
therefore for pure compatibility -P and -Q should rather be sparse
than full (and zeros should be +0 rather than -0) and no one would
be really surprised if Octave silently converts the permutation
matrices to sparse, whereas some code may be negatively impacted
if they are converted to full and generate a memory overflow.

Carlo de Falco <cdf> Fri 25 Mar 2016 02:19:20 AM UTC, comment #25:

It was Carlo, not me, who suggested the sign bit.

I don't see the need to get rid of permutation and diagonal matrices over this issue.  (As you say, code size is another matter entirely.)

I'm just worried about sequences like

A = eye(3);
A = -A(end:-1:1, :);
A(1,1,2) = 1;

My understanding is that if the second command causes A to be sparse, then the third command will cause an error.  That level of incompatibility is problematic.  Getting the sign of zero to agree with Matlab is much less important, even if people then divide by it.

Is it documented anywhere how diagonal and permutation matrices are formed?  I assume that diagonal matrices are much more common.  (Contrary to what Carlo said, permutations are not a generalisation of diagonal matrices, since permutations are 0/1 matrices.) Rather than adding a sign to a permutation matrix, how about allowing exactly one arbitrary value in each row/column.  That way they are actually generalisations of diagonal matrices.  They would be slower and consume more memory than now, but that may be more than compensated for by being more widely used.

Another option would be to drop permutation matrices, and keep diagonal matrices.  That presumably saves half the lines of code, but still optimises a very common class of matrices.  Out of interest, why were permutation matrices chosen to optimise, rather than, say, tridiagonal matrices, triangular matrices or Hermitian matrices?  They all seem more common and allow more efficient algorithms to be chosen.

Lachlan Andrew <lachlan>
Thu 24 Mar 2016 06:21:51 PM UTC, comment #24:

I don't see an easy solution to perfect Matlab compatibility and also having Octave language extensions.  Actually, there is an easy solution which is to disable extensions when perfect compatibility is required, but jwe wasn't keen on that.

If we require exact compatibility, then we might want to just get rid of both Permutation and Diagonal matrices entirely.  They do add a lot of code in terms of files, operators between P & D matrices and every other type, and in switch statements where we have to decode the matrix type.

But we would be losing a lot in time and memory.  For example,

octave:28> x = rand (1e3);
octave:29> P = eye (1e3)(1000:-1:1, :);
octave:30> Ps = sparse (P);
octave:31> Pf = full (P);
octave:32> tic; y = P*x; toc
Elapsed time is 0.00626493 seconds.
octave:33> tic; ys = Ps*x; toc
Elapsed time is 0.0852761 seconds.
octave:34> tic; yf = Pf*x; toc
Elapsed time is 0.114365 seconds.
octave:35> whos P*
Variables in the current scope:

Attr Name        Size                     Bytes  Class
==== ====        ====                     =====  =====
P        1000x1000                    4000  double
Pf       1000x1000                 8000000  double
Ps       1000x1000                   16004  double

The Permutation matrix is the smallest by far, and executes 13X faster then the Sparse matrix implementation.  The full matrix implementation is not only tremendously wasteful of space, it also is another 30% slower than the Sparse implementation.

I suppose if it was super important to preserve '-0' then we could add a signbit flag to Permutation and Diagonal matrices to record that fact as suggested by Lachlan.

Rik <rik5> Thu 24 Mar 2016 10:33:24 AM UTC, comment #23:

Lachlan,

Maybe I missed something, but I cannot think of a case when a user creates a full matrix and Octave converts it to permutation matrix automatically, did you ever see this happen? What is an example?

I was assuming that automatic conversions only occur from the most specialized to the most general class, i.e.

diagonal -> permutation -> sparse -> full

but never the other way, is this assumption wrong?

Carlo de Falco <cdf> Thu 24 Mar 2016 05:43:01 AM UTC, comment #22:

That is true, Carlo, but it introduces a clear incompatibility with Matlab.

If a user creates a full matrix, and Octave converts that to a permutation matrix, and then to a sparse matrix, what happens when the user tries an operation that is valid only for full matrices (like assigning to a third dimension)?

In Matlab, the operation would succeed because the matrices would always be full.

Lachlan Andrew <lachlan>
Thu 24 Mar 2016 05:31:55 AM UTC, comment #21:

JWE I think exposing the details of the permutation matrix class
implementation to the end user would be really a bad thing.

The purpose of using a special class for permutation matrices is to exploit the benefits in terms of speed and memory efficiency that are available in this case.

For the case of a row/column permuted diagonal matrix, the next best possible implementation in terms of speed and memory efficiency is that of a sparse matrix.

Therefore converting a permutation matrix to sparse whenever an operation that is not implemented is requested makes much more sense than falling back to full matrices.

I agree with rik that low level details like handling overflows and '-0' should be given lower priority.

Carlo de Falco <cdf> Wed 23 Mar 2016 08:59:11 PM UTC, comment #20:

How about limiting the operations that are available with  permutation matrices matrices so that they can only permute the rows or columns of other matrices?  Then this issue disappears and you will have to explicitly convert to sparse or full if you want to perform other operations on them.  Or, disable the special feature.

John W. Eaton <jwe> Wed 23 Mar 2016 08:02:41 PM UTC, comment #19:

The numerical error is large only if you define it in reference to Matlab, rather than in reference to pure math.  In mathematics, 0 is unsigned so 1 / -0 == Inf.

It seems like we are allowing the details of the implementation (the fact that the IEEE floating point standard has a signbit separate from the actual number) to dictate what the result should be rather than looking at the specification (the rules of mathematics).  Normally, leaking of the implementation through an API is considered a bad thing as the results may change wildly based on the underlying implementation.

I'm going to use Diagonal matrices for examples which have the same issues as Permutation matrices.  Diagonal, Permutation, and Sparse matrices will all need to be reviewed after a path forward is chosen.

x = diag ([1 2])
x =

Diagonal Matrix

1   0
0   2

That's okay, but now take the negation

y = -x
y =

Diagonal Matrix

-1   0
0  -2

This agrees with mathematics and with sparse matrices, but disagrees with Permutation matrices and Full matrices.

For regular Full matrices,

xd = [1 0; 0 2]
xd =

1   0
0   2

yd = -x
yd =

-1  -0
-0  -2

But this is implementation dependent.  If the numbers are stored as two's complement then there is only one zero and the behavior, even for Full matrices, is inconsistent.

xi = int8 ([1 0; 0 2])
xi =

1  0
0  2

yi = -xi
yi =

-1   0
0  -2

In quantum mechanics you often get answers that are the same up to an arbitrary phase factor.  It seems like this is the case here.  The value is always zero, but the phase factor is +/- 1.  With multiplication, the value is 0 so the result will be be +/- 0.  With division, the value is zero so the result will be +/- Inf.

I don't think it is unreasonable to have Octave pick the arbitrary phase factor in a way that it easy for us to understand (just like pure math) or to code.  If programmers want exact congruence with Matlab then they can pay for Matlab or disable Octave's special features.

Rik <rik5> Wed 23 Mar 2016 05:24:43 PM UTC, comment #18:

I don't really like the idea of getting different numerical results based on configuration variables.  I mean, it would be OK to be off by a few eps, but -Inf vs +Inf is about as big a numerical error as you can get.

John W. Eaton <jwe> Wed 23 Mar 2016 05:13:51 PM UTC, comment #17:

It seems like maybe we are worrying overly much about corner cases.  What about implementing things in Octave in a space efficient manner, and if people absolutely need Matlab compatibility then have them use

disable_permutation_matrix (true)

Rik <rik5> Wed 23 Mar 2016 04:59:23 PM UTC, comment #16:

Yes, Matlab does not have permutation matrices.  So I think the question is, do we make the result of things like

a = rand (2);
[q, r, p] = qr (a);
neg_p = -p;
1 ./ neg_p;

produce the same result as Matlab, or not?

Even though Matlab will not display the -0, I think it will produce -Inf for the division by -0 because it does compute the -0, it's just not displaying it directly.

John W. Eaton <jwe> Wed 23 Mar 2016 04:48:35 AM UTC, comment #15:

Matlab doesn't have Diagonal matrices in the same way that Octave does.  Octave has special code that only stores the elements on the diagonal, tremendously reducing storage requirements; While Matlab just implements things as a full matrix.

As for the question of equality,

-full (P) == full (-P)

this relationship will hold even if Permutation matrices do not implement -0 because

octave:1> 0 == -0
ans =  1

Rik <rik5> Wed 23 Mar 2016 12:58:09 AM UTC, comment #14:

On 8.2.0.701 (R2013b), it produces

testsignbit (1)
testsignbit (-1)
testsignbit (0)
testsignbit (-0)
Hello, World!
Hello, World!
SIGNBIT!
Hello, World!
Hello, World!
SIGNBIT!

It seems that this extends to retaining the sign of zeros in diagonal matrices:

>> d = diag ([1, 2])
d =
1     0
0     2
>> -d
ans =
-1     0
0    -2
>> 1./-d
ans =
-1.0000      -Inf
-Inf   -0.5000

Comment #5 asked if it is important that -full (P) and full (-P) produce the same result.  I think it is important that

-full (P) == full (-P)

but don't think it is important the result be identical.

However, if permutation matrices are only automatically created from full matrices, then it is probably not too bad if they are implicitly converted back.  If the sparsity is required, it would be easy to type

-sparse(jacQ)

Does Octave ever create a permutation matrix from a sparse matrix?

Lachlan Andrew <lachlan>
Tue 22 Mar 2016 09:49:09 PM UTC, comment #13:

Are you sure that it throws away the sign on -0?  Maybe it is just not displayed?

Could someone try the attached MEX file with Matlab and see what it does for the following calls?

testsignbit (1)
testsignbit (-1)
testsignbit (0)
testsignbit (-0)

John W. Eaton <jwe> Tue 22 Mar 2016 08:58:32 PM UTC, comment #12:

Answering an earlier question: Matlab does not distinguish between 0 and -0 and stores everything as +0.  Octave, for full matrices, does have both +0 and -0.  Sparse matrices, however, only store +0.

In my view, sparse, permutation, and diagonal matrices are all closer to each other than they are to full matrices.  Ergo, I also believe that permutation matrices should behave more like sparse matrices under negation.  In this case, -P should only negate the non-zero elements of P.  If we're looking for consistency, this is already what Diagonal matrices do.

z = diag ([1 2 3])
z =

Diagonal Matrix

1   0   0
0   2   0
0   0   3

octave:11> z2 = -z
z2 =

Diagonal Matrix

-1   0   0
0  -2   0
0   0  -3

Rik <rik5> Mon 21 Mar 2016 10:24:27 AM UTC, comment #11:

I was considering whether it would make sense
to add a 'sign' property to the permutation matrix
class in order to avoid the conversion altogether.

But after running the following test:

n = 1e4;
a = eye (n);
a = a ([end 2:end-1 1], :);
typeinfo (a)
b = sparse (a);
typeinfo (b)
c = sprandn (n);
tic, a * c; toc
tic, b * c; toc
c = randn (n);
tic, b * c; toc
tic, a * c; toc

I see that matrix-matrix multiplication is
consistently at least twice as fast with sparse
matrices than permutation matrices, so I
think lachlan's idea makes more sense ...

c.

Carlo de Falco <cdf> Mon 21 Mar 2016 01:14:07 AM UTC, comment #10:

I've changed the title (from "out of memory with factorized sparse matrix using mldivide") to reflect jwe's explanation.

IMHO, negating a permutation matrix should give a sparse matrix.  Distinguishing between 0 and -0 is not mathematically sound in the first pace, and not worth forcing a sparse matrix to be full.

Lachlan Andrew <lachlan>
Tue 15 Mar 2016 09:26:26 PM UTC, comment #9:

My last two posts don't show up correctly in the tracker
so I'm trying to post them again.

Octave 4.1

[l, u, p] = lu (A);
issparse (p)
ans = 0
issparse (-p)
ans = 0
-verbaitim-

Matlab

+verbatim+
[l, u, p] = lu (A);
issparse (-p)
ans =
1
issparse (p)
ans =
1
-verbaitim-

In Octave 3.8.1 I get:

+verbatim+
octave:4> A = sprandn (10, 10, .4);
octave:5> [l, u, p] = lu (A);
octave:6> issparse (p)
ans = 0
octave:7> issparse (-p)
ans =  1
octave:8> version
ans = 3.8.1

this latter seems the most correct to me ...
is there a reason why it was changed?

Carlo de Falco <cdf> Tue 15 Mar 2016 09:18:35 PM UTC, comment #8:

Octave 4.1+

>> [l, u, p] = lu (A);
>> issparse (p)
ans = 0
>> issparse (-p)
ans = 0
-verbaitim-

Matlab

+verbatim+
>> [l, u, p] = lu (A);
>> issparse (-p)
ans =
1
>> issparse (p)
ans =
1
-verbaitim-

In Octave 3.8.1 I get:

+verbatim+
octave:4> A = sprandn (10, 10, .4);
octave:5> [l, u, p] = lu (A);
octave:6> issparse (p)
ans = 0
octave:7> issparse (-p)
ans =  1
octave:8> version
ans = 3.8.1

this latter seems the most correct to me ...
is there a reason why it was changed?

Carlo de Falco <cdf> Tue 15 Mar 2016 08:59:31 PM UTC, comment #7:

I also tried the following:

Octave 4.1+

>> [l, u, p] = lu (A);
>> issparse (p)
ans = 0
>> issparse (-p)
ans = 0
-verbaitim-

Matlab
+verbatim+
>> [l, u, p] = lu (A);
>> issparse (-p)
ans =
1
>> issparse (p)
ans =
1
-verbaitim-

this difference makes sense because AFAIK Matlab, unlike Octave
does not have a specific permutation matrix class, but should we also return true for issparse when applied to a permutation (or diagonal) matrix? should unary minus applied to a permutation matrix return a sparse rahter than full matrix?

Carlo de Falco <cdf> Tue 15 Mar 2016 08:52:00 PM UTC, comment #6:

I don't know whether it is relevant, but I tried the following in Octave and Matlab:

>> version
ans =
8.6.0.267246 (R2015b)
>> A = sprandn (5, 5, .4);
>> full (A)
ans =
0.3035         0         0   -0.6003         0
0.4900         0         0         0         0
0    0.7394         0         0         0
0         0    1.7119         0         0
-0.1941         0         0         0         0
>> full (-A)
ans =
-0.3035         0         0    0.6003         0
-0.4900         0         0         0         0
0   -0.7394         0         0         0
0         0   -1.7119         0         0
0.1941         0         0         0         0

>> -full (A)
ans =
-0.3035         0         0    0.6003         0
-0.4900         0         0         0         0
0   -0.7394         0         0         0
0         0   -1.7119         0         0
0.1941         0         0         0         0

>> version
ans = 4.1.0+
>> A = sprandn (5, 5, .4);
>> full (A)
ans =
-0.74912   0.00000   0.00000   0.00000   0.00000
0.00000   0.00000   0.00000   0.00000   0.00000
0.00000   0.00000  -0.54435   1.02264   0.00000
0.00000  -0.95284   0.00000  -0.64702   0.21304
-0.78968   0.49972  -1.67818   0.00000  -0.07748
>> full (-A)
ans =
0.74912   0.00000   0.00000   0.00000   0.00000
0.00000   0.00000   0.00000   0.00000   0.00000
0.00000   0.00000   0.54435  -1.02264   0.00000
0.00000   0.95284   0.00000   0.64702  -0.21304
0.78968  -0.49972   1.67818   0.00000   0.07748
>> -full (A)
ans =
0.74912  -0.00000  -0.00000  -0.00000  -0.00000
-0.00000  -0.00000  -0.00000  -0.00000  -0.00000
-0.00000  -0.00000   0.54435  -1.02264  -0.00000
-0.00000   0.95284  -0.00000   0.64702  -0.21304
0.78968  -0.49972   1.67818  -0.00000   0.07748
>>

does this mean Matlab returns 0 for negating the zero entries in a full matrix while Octave returns -0 ?

Carlo de Falco <cdf> Tue 15 Mar 2016 06:52:42 PM UTC, comment #5:

The problem isn't the multiplication or division steps, it is the use of the unary minus operator on jacQ.  Your single line expression isn't exactly the same as the series of operations.  To be equivalent, you need to write

delta = - (jacQ * (jacU \ (jacL \ (jacP * (jacR \ res)))));

You will also see the same error if you write

-jacQ

The reason is that jacQ is a permutation matrix, not a sparse matrix, and for permutation matrices, negation returns -0 for all the 0 elements of the original permutation matrix and that forces it to be converted to a full matrix.  For a sparse matrix, the unary minus operator ignores the 0 elements.

I'm not sure what the best solution is.  I suppose the real question is whether the following should produce the same result for a permutation matrix P:

-full (P)
full (-P)

Currently this is true for permutation matrices, but not for sparse matrices.

John W. Eaton <jwe> Tue 15 Mar 2016 12:43:37 PM UTC, comment #4:

sorry I had a typo in the previous post the working code is:

>> a = jacR \ res;
>> b = jacP * a;
>> c = jacL \ b;
>> d = jacU \ c;
>> e = jacQ * d;
>> delta = -e;

while doing the same on one line results in an error

>> delta = - jacQ * (jacU \ (jacL \ (jacP * (jacR \ res))));
error: out of memory or dimension too large for Octave's index type

Carlo de Falco <cdf> Tue 15 Mar 2016 10:47:15 AM UTC, comment #3:

Intrestingly enough, I tried separating the solve line into
different steps and it works:

>> [jacL, jacU, jacP, jacQ, jacR] = lu (jac);
>> delta = - jacQ * (jacU \ (jacL \ (jacP * (jacR \ res))));
error: out of memory or dimension too large for Octave's index type
>> a = jacR \ res;
>> b = jacP \ a;
>> c = jacL \ b;
>> d = jacU \ c;
>> e = jacQ * d;
>> delta = -e;

maybe the error lies in the parser somewhere?

Carlo de Falco <cdf> Tue 15 Mar 2016 10:42:08 AM UTC, comment #2:

probably the file is too big for the tracker.
https://www.dropbox.com/s/jqligynj8ffoy91/sparse_matrix_bug_data.binary.gz?dl=0

c.

Carlo de Falco <cdf> Tue 15 Mar 2016 09:39:11 AM UTC, comment #1:

Carlo, I can't see the attachment.  Could you please try again?

You ask for suggestions of how to debug it.  The natural thing would be to work out which particular operation in the last line is failing.

Lachlan Andrew <lachlan>
Tue 15 Mar 2016 09:23:46 AM UTC, original submission:

Hi,

The following test with a moderately large
sparse matrix results in an error on the release
branch but it works on version 3.8.1, unfortunately
I do not have the stable version at hand at the moment,
can someone check what happens there?

>>[jacL, jacU, jacP, jacQ, jacR] = lu (jac);
>>delta = - jacQ * (jacU \ (jacL \ (jacP * (jacR \ res))));
error: out of memory or dimension too large for Octave's index type

the matrix that caused this is attached.
Can anyone suggest how to further debug this?

c.

Carlo de Falco <cdf> ## Attached Files

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

Attach Files:

Comment:

Attached Files
file #36727:  testsignbit.c added by jwe (291B - text/x-csrc)

## Dependencies

Depends on the following items: None found

Items that depend on this one: None found

Carbon-Copy List
• -email is unavailable- added by nrjank (Posted a comment)
• -email is unavailable- added by rik5 (Posted a comment)
• -email is unavailable- added by jwe (Posted a comment)
• -email is unavailable- added by lachlan (Posted a comment)
• -email is unavailable- added by cdf (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 project members can vote.

## History

2016-03-22 jwe Attached File- Added testsignbit.c, #36727
2016-03-21 lachlan Summaryout of memory with factorized sparse matrix using mldivide out of memory negating a permutation matrix