Re: gives different result compared to 1/Diagonal[Normal@A] when A is sparse
- To: mathgroup at smc.vnet.net
- Subject: [mg123342] Re: gives different result compared to 1/Diagonal[Normal@A] when A is sparse
- From: Leonid Shifrin <lshifr at gmail.com>
- Date: Fri, 2 Dec 2011 07:21:23 -0500 (EST)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
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 <danl at 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 > >> <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 >>>> >>>> >>>> >>> >>> >