Re: gives different result compared to 1/Diagonal[Normal@A] when A is sparse

[Oliver Ruebenkoenig]

I filed this as a bug.

One small side note: You might want to use

A[[r, r + 1]] = 0.;
A[[r + 1, r]] = 0.;

instead of A[[..]]=0.& /@ r. Extraction/Setting of SA componets is
vectoriezed.

Oliver

On Wed, 30 Nov 2011, Nasser M. Abbasi wrote:

> I found this strange behavior, and I do not think it is correct.
>
> This is version 8.04.
>
> 1/Diagonal[A] gives a divide by zero error, but 1/Diagonal[Normal@A] does
> not. This is when A is sparse.
>
> ------------------------------
> Clear["Global`*"]
>
> makeMatrix[n_]:=Module[{numberOfUnknowns=n^2,r,A},
>
> A=SparseArray[
> {
> Band[{1,1}]->4.0,
> Band[{2,1}]->-1,
> Band[{1,2}]->-1,
> Band[{1,n+1}]->-1,
> Band[{n+1,1}]->-1
> },{numberOfUnknowns,numberOfUnknowns},0.
> ];
>
> r=Range[n,n^2-n,n];
> (A[[#,#+1]]=0.)&/@r;
> (A[[#+1,#]]=0.)&/@r;
>
> A
> ];
>
> (A = makeMatrix[3])//MatrixForm
>
> (Diagonal[A])//Normal
>
> 1/Diagonal[Normal@A] (* ===> OK *)
> 1/Diagonal[A] (* error *)
>
> ----------------------------------------
>
> So, 1/Diagonal[Normal@A] gives
> {0.25,0.25,0.25,0.25,0.25,0.25,0.25,0.25,0.25}
>
> but 1/Diagonal[A] gives 1/0
>
> In another system I use, both operations give
> {0.25,0.25,0.25,0.25,0.25,0.25,0.25,0.25,0.25}
>
> i.e. if the matrix is sparse or not, 1/Diagonal[A] should work
> regardless. I think sparse matrices need to be more integrated into
> all Mathemaitca matrix operations.
>
> Or Am I missing something here?
>
> Thanks,
> --Nasser
>
>

[Leonid Shifrin]

I think this is the same bug as discussed in this thread:

http://groups.google.com/group/comp.soft-sys.math.mathematica/browse_thread/thread/37935f55b291ed2

where in my answer there I also indicated what the origin of that behavior
might be.

Cheers,
Leonid

On Wed, Nov 30, 2011 at 3:50 PM, Oliver Ruebenkoenig
<rueb...@wolfram.com>wrote:

[Leonid Shifrin]

Daniel,

Upon some reflection, I think I agree with you. The only thing I'd improve is
the wording of the error messages. Perhaps some custom error message
could be issued in such cases - currently the message gives little clue for
what could be causing it. This looks almost like a leak of implementation
detail into the public interface. I think, SparseArray is under-specified
in
this regard. It should be stated in the docs that when a reciprocal is taken
of the SparseArray, the default element is the inverse of the original
default
element, including cases when the latter is zero. And the error message
should be different in this case, giving a specific warning about default
element, rather than a general 1/0 warning which looks like a mystery,
especially in cases like

SparseArray[{1}]/SparseArray[{1}]

where no explicit zeros are involved whatsoever.


Cheers,
Leonid


On Thu, Dec 1, 2011 at 7:31 PM, Daniel Lichtblau <da...@wolfram.com> wrote:

> On 12/01/2011 04:49 AM, Leonid Shifrin wrote:
>
>> I think this is the same bug as discussed in this thread:
>>

>> http://groups.google.com/**group/comp.soft-sys.math.**
>> mathematica/browse_thread/**thread/37935f55b291ed2<http://groups.google.com/group/comp.soft-sys.math.mathematica/browse_thread/thread/37935f55b291ed2>

>>
>> where in my answer there I also indicated what the origin of that behavior
>> might be.
>>
>> Cheers,
>> Leonid
>>
>

> I think the current behavior is The Right Thing, in a mathematical sense
> (not due to structural reasons such as Listable attribute). If one does not
> invert the value of implied element, then what would happen in the example
> below? (We start with the definitions from the original post, found below.)
>
> aa = makeMatrix[3];
>
> In[57]:= bb = 1/Diagonal[aa];
>
> During evaluation of In[57]:= Power::infy: Infinite expression 1/0.
> encountered. >>
>
> In[58]:= cc = SparseArray[ArrayRules[bb] /. {xx_, yy___} -> {yy}];
>
> In[59]:= Normal[cc]
> Out[59]=
> {ComplexInfinity, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25}
>
> Compare to reversing the element removal with reciprocation.
>
> In[60]:= aa2 = Diagonal[aa];
>
> In[61]:= bb2 = SparseArray[ArrayRules[aa2] /. {xx_, yy___} -> {yy}];
>
> In[62]:= cc2 = 1/bb2;
>
> During evaluation of In[62]:= Power::infy: Infinite expression 1/0.
> encountered. >>
>
> In[63]:= Normal[cc2]
>
> Out[63]= {ComplexInfinity, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, \
> 0.25}
>
> I think it is good that these agree. But that would not happen, were there
> a special case to avoid inverting the implied element when all possible
> elements are explicitly present. Stated more generally, such a special case
> would give rise to inconsistancies.
>
> Possibly the Power::infy message should be suppressed. That would be a
> minor bug, if it is a bug at all (more of a design decision, I think).
>
> Daniel Lichtblau
> Wolfram Research

>
>
>
> 0, 2011 at 3:50 PM, Oliver Ruebenkoenig
>
>> <rueb...@wolfram.com>wrote:
>>
>>
>>>
>>> I filed this as a bug.
>>>
>>> One small side note: You might want to use
>>>
>>> A[[r, r + 1]] = 0.;
>>> A[[r + 1, r]] = 0.;
>>>
>>> instead of A[[..]]=0.& /@ r. Extraction/Setting of SA componets is
>>> vectoriezed.
>>>
>>> Oliver
>>>
>>> On Wed, 30 Nov 2011, Nasser M. Abbasi wrote:
>>>
>>> I found this strange behavior, and I do not think it is correct.
>>>>
>>>> This is version 8.04.
>>>>
>>>> 1/Diagonal[A] gives a divide by zero error, but 1/Diagonal[Normal@A]
>>>>
>>> does
>>>
>>>> not. This is when A is sparse.
>>>>
>>>> ------------------------------
>>>> Clear["Global`*"]
>>>>

>>>> makeMatrix[n_]:=Module[{**numberOfUnknowns=n^2,r,A},

>>>>
>>>> A=SparseArray[
>>>> {
>>>> Band[{1,1}]->4.0,
>>>> Band[{2,1}]->-1,
>>>> Band[{1,2}]->-1,
>>>> Band[{1,n+1}]->-1,
>>>> Band[{n+1,1}]->-1

>>>> },{numberOfUnknowns,**numberOfUnknowns},0.

>>>> ];
>>>>
>>>> r=Range[n,n^2-n,n];
>>>> (A[[#,#+1]]=0.)&/@r;
>>>> (A[[#+1,#]]=0.)&/@r;
>>>>
>>>> A
>>>> ];
>>>>
>>>> (A = makeMatrix[3])//MatrixForm
>>>>
>>>> (Diagonal[A])//Normal
>>>>
>>>> 1/Diagonal[Normal@A] (* ===> OK *)
>>>> 1/Diagonal[A] (* error *)
>>>>

>>>> ------------------------------**----------
>>>>
>>>> So, 1/Diagonal[Normal@A] gives
>>>> {0.25,0.25,0.25,0.25,0.25,0.**25,0.25,0.25,0.25}

>>>>
>>>> but 1/Diagonal[A] gives 1/0
>>>>
>>>> In another system I use, both operations give

>>>> {0.25,0.25,0.25,0.25,0.25,0.**25,0.25,0.25,0.25}