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