MathGroup Archive 2011

[Date Index] [Thread Index] [Author Index]

Search the Archive

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

  • To: mathgroup at smc.vnet.net
  • Subject: [mg123351] Re: gives different result compared to 1/Diagonal[Normal@A] when A is sparse
  • From: Daniel Lichtblau <danl at wolfram.com>
  • Date: Fri, 2 Dec 2011 07:23:01 -0500 (EST)
  • Delivered-to: l-mathgroup@mail-archive0.wolfram.com

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
>
> 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
> <ruebenko at 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}
>>>
>>> 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
>>>
>>>
>>
>>




  • Prev by Date: Re: Problem with Coding for Box Whisker Chart
  • Next by Date: Re: Orthogonalize[expr]
  • Previous by thread: Re: gives different result compared to 1/Diagonal[Normal@A] when A is sparse
  • Next by thread: Re: A function to do incomplete LU decomposition with a drop tolerance?