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?

issparse(P) returns zero.
in the workspace, P is listed as a 100x100 double.

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

thought I'd add some info but both :
produce the error

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,
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
Also, if I understand things correctly it shouldn't return the Q and R values as those are only for sparse matrices.

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:
In Octave 4.1 the same code returns
while the conversion to sparse would lead to the following
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??

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:
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
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.

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
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.

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,
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.

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?

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.

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.

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.

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.
That's okay, but now take the negation
This agrees with mathematics and with sparse matrices, but disagrees with Permutation matrices and Full matrices.
For regular Full matrices,
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.
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.

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.

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

Yes, Matlab does not have permutation matrices. So I think the question is, do we make the result of things like
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.

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,
this relationship will hold even if Permutation matrices do not implement 0 because

On 8.2.0.701 (R2013b), it produces
It seems that this extends to retaining the sign of zeros in diagonal matrices:
Comment #5 asked if it is important that full (P) and full (P) produce the same result. I think it is important that
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
Does Octave ever create a permutation matrix from a sparse matrix?

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?
(file #36727)

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 nonzero elements of P. If we're looking for consistency, this is already what Diagonal matrices do.

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:
I see that matrixmatrix multiplication is
consistently at least twice as fast with sparse
matrices than permutation matrices, so I
think lachlan's idea makes more sense ...
c.

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.

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

Octave 4.1+

I also tried the following:
Octave 4.1+

I don't know whether it is relevant, but I tried the following in Octave and Matlab:
does this mean Matlab returns 0 for negating the zero entries in a full matrix while Octave returns 0 ?

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
You will also see the same error if you write
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:
Currently this is true for permutation matrices, but not for sparse matrices.

sorry I had a typo in the previous post the working code is:
while doing the same on one line results in an error

Intrestingly enough, I tried separating the solve line into
different steps and it works:
maybe the error lies in the parser somewhere?

probably the file is too big for the tracker.
You can download it from here:
https://www.dropbox.com/s/jqligynj8ffoy91/sparse_matrix_bug_data.binary.gz?dl=0
c.

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.

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?
the matrix that caused this is attached.
Can anyone suggest how to further debug this?
c.
