Wed 07 Oct 2015 02:46:26 PM UTC, comment #17:
I changed the behavior of Octave to return Inf where the diagonal element is 0 and to preserve the return class as diagonal matrix. See this changeset (http://hg.savannah.gnu.org/hgweb/octave/rev/f61c67865d9f). To access the fix the original reporter will need to build from source, or wait until release 4.2.
|
Sat 03 Oct 2015 01:46:27 PM UTC, comment #16:
Hi,
I teach my students that the inverse of a diagonal matrix is simply a diagonal matrix with the inverse diagonal. So, I prefer solution a) in #6, I consider it the most natural and the least surprising. I agree with #5 that Octave should give the same result as with a regular matrix, but the challenge should be have Octave to recognize a regular diagonal matrix A as a diagonal matrix and compute the inverse by diag(1./diag(A)). My 2 cents...
Marco
|
Sat 03 Oct 2015 09:09:29 AM UTC, comment #15:
A sidenote on the questions of sparsity:
blkdiag() returns sparse if one of the arguments is sparse. I believe diagonal matrices - at least here - should be considered sparse rather that full? Is there a general rule on which side diagonal matrix falls? I wonder if issparse() should return true for diagonal matrix?
|
Sat 03 Oct 2015 06:42:08 AM UTC, comment #14:
Thank you everybody for quick replies. Still some comments:
Indeed, diagonal (and e.g. permutation) matrices are one of the few, but important, classes that (mathematically) preserve the very sparsity they posess under inversion. In fact diagonal matrices form a group, so all algebraic operations preserve this structure. I think this should be respected in the code.
Now, I, too, tried this in matlab with the following finndings
1) diag() produces a full matrix, i.e. no respective diagonal class exists.
2) Under the inversion of the sparse matrix a sparse result is produced.
3) the inversion of a spoarse matrix seems to be some Gaussian type of elimimination with potential divisions by zero producing inf and divisions by inf nan, without any warnings. This makes sense and I don't think this is a bad solution.
I would be more forgiving to matlab. As a matter of principle: Producing infs and nans as a natural product of (elementwise) operations (when feasible) follows the principle of "least loss of information" and indeed makes the debugging of operations in the code easiest. So I would still keep my vote for alternative (a) in post #6.
|
Fri 02 Oct 2015 11:46:27 PM UTC, comment #13:
WTF Matlab? That is a weird result for the sparse case. I don't think we should follow that. I'll look at modifying inv.cc.
|
Fri 02 Oct 2015 11:34:40 PM UTC, comment #12:
Some Matlab 2015a results:
|
Fri 02 Oct 2015 11:25:53 PM UTC, comment #11:
I'll ask on the Maintainer's list for someone to run the Matlab code. I knew the mutation would be hard so I think we should avoid it in dDiagMatrix.cc. Maybe, as I said, we can figure it out in inv.cc and do the mutation there.
But, I don't think we should always promote to full. Lots, even most, of the diagonal matrices will not be singular and so we don't want to go from linear memory consumption (N) to quadratic consumption (N^2) every time.
|
Fri 02 Oct 2015 11:18:45 PM UTC, comment #10:
Sorry, that's it for my Matlab hookup for now, someone else will have to test the sparse case.
I started working on a changeset for this, but changing the functions so that they return either a diagonal matrix or a full matrix proved a bit more involved than I thought.
If the inverse of a diagonal always returns a full matrix, singular or not, do we think that is also an acceptable solution?
|
Fri 02 Oct 2015 09:56:50 PM UTC, comment #9:
@Mike: Could you try the sparse example under Matlab?
Under Octave it simply refuses to calculate an answer
This could be another option to use for singular diag matrices. It would certainly be easier to implement. The C++ inverse function defined in dDiagMatrix.cc returns a DiagMatrix, not an octave_value. Thus we can't easily mutate the return value from the special DiagMatrix class to the ordinary Matrix class. I suppose one could do so at a higher level in inv.cc by checking for any element equal to 0. That's pretty easy check on a diagonal matrix.
|
Fri 02 Oct 2015 09:14:29 PM UTC, comment #8:
Ok, the argument of compatibility (with matlab) is valid. Does matlab nowadays implement diagonal class? If not matlab's behaviour might change in future.
What is behaviour of matlab for the case of the sparse matrixes.
Apart from the current behaviour, all mentioned solutions are acceptable for a mathematician (at least this one). Still the case (a) in post #6 however is the most natural one to me.
|
Fri 02 Oct 2015 09:10:44 PM UTC, comment #7:
Indeed in Matlab it returns a full matrix of inf. Despite the possibility of memory overflow if the user is expecting the operation to return a diag of the same size as the input, I think that's the best fix to make here.
|
Fri 02 Oct 2015 09:07:09 PM UTC, comment #6:
a correction to my previous post:
My suggestion is either a) or b) below, due to the potential memory blow up.
A=diag([1,0])
a) inv(A) -> diag([1,inf])
b) inv(A) -> diag(inf,inf])
c) inv(A) -> [inf,inf;inf,inf];
|
Fri 02 Oct 2015 09:04:01 PM UTC, comment #5:
For something simple like 2+2 there is a very well-defined output. For something more complex like the inverse of a singular matrix there are several possibilities and arguments can be made for each one. I am advocating following the Principle of Least Surprise. The implementation of diagonal matrices is optimized in Octave, but users don't need to know or understand that. To them, a diagonal matrix is the same as a regular matrix and the same calculation should yield the same result irrespective of what Octave has used for the underlying object type.
One useful test is to have the following code run under Matlab
If this produces a 3x3 matrix of Inf values then that is what Octave should also produce.
I can't tell from comment #3 which is considered the odd behavior: 1) returning the original matrix, or 2) returning an all Inf matrix. I certainly want to stop using method 1. And for the next question, I would prefer to return Inf (3,3) rather than diag (Inf (3, 1))
|
Fri 02 Oct 2015 08:59:00 PM UTC, comment #4:
Indeed, there is already difference in the behaviour of the inverse operation for the sparse and respectively the full matrix classes. The sparse matrix only gives an error and no result as shown below.
One aspect of using the elementwise inversion (with 1/0=inf) is that it preserves the sparsity structure and does not produce memory blow up. If one wants to follow the approach similar to the full matrix case, perhaps a diagonal matrix with the entire column inf's is a better choice:
inv(diag([1,0])) -> diag([1,inf])
|
Fri 02 Oct 2015 08:39:55 PM UTC, comment #3:
With all due respect, inverse is mathematically a different concept to pseudo-inverse. And no function should knowlingly give such a misleading result that can give very odd behaviour.
And since the diagonal matrix is now produced by e.g. diag(), already the bw-compatibility requires this.
It seems that simply the following patch would do, but as this is my first peak into octave source, I might be wrong.
|
Fri 02 Oct 2015 03:49:52 PM UTC, comment #2:
The diagonal matrix class is supposed to be a memory saving optimization. All calculations on diagonal matrices should result in exactly the same behavior as when a full matrix is used. My suggestion would be to change the diagonal class to return an Inf array for a singular matrix the same as is done fore a full matrix.
|
Fri 02 Oct 2015 01:12:58 PM UTC, comment #1:
Thanks for your bug report. This does not appear to be documented, but when inverting a singular diagonal matrix, the original matrix is returned.
As you've noticed, and as the documentation for inv() clearly states, there are often better ways to solve the problem using mldivide, or pinv.
Also, if you like this answer better,
|
Fri 02 Oct 2015 09:51:36 AM UTC, original submission:
Inverting (with inv) a singular diagonal matrix results the
matrix itself, but mrdivide and mldivide act correctly.
This bug affects also at least version 3.8.2. (from ubuntu 15.04. packaging.)
|