- petsc-dev - mailman3.cels.anl.gov

Re: ufuncs, iterators
by Barry Smith 19 Aug '05

19 Aug '05

After sleeping on it, it may be ok to include these methods in Mat as methods for "filling up matrices" so long as in the end you end up with a Mat that you then use as an operator. But I'd still like to see/understand a little more of the "construction" process. Classically one would do that as loops over elements and perform all the computations for the one element before moving to the next. In Matlab this can be done instead (with some impact on performance) using a sequence of array operations (with the loop inside each array operation). In the past, since PETSc was used exclusively from C/C++ and Fortran users built their matrices (operators) using the "classical" approach, now with python it appears reasonable that we may need to add the "array" approach. Barry On Fri, 19 Aug 2005, Matthew Knepley wrote: > Simon Burton <simon(a)arrowtheory.com> writes: > > > On Thu, 18 Aug 2005 19:57:02 -0500 > > Matthew Knepley <knepley(a)mcs.anl.gov> wrote: > > > >> > >> > (d) m=exp(m) (pointwise exp) > >> > >> We could add pointwise operations just like the VecPointwise*(). > >> > >> Matt > > > > Yes, we need MatPointwiseMult aswell. > > > > Should I make a start on this ? > > I am thinking more about what Barry said. The VecPointwise*() operations > can be given a solid mathematical interpretation in terms of spinor operations. > However, I do not see anything like that for the Mat stuff yet. We need to > understand the mathematicas better. > > Matt >

1 0

Re: ufuncs, iterators
by Matthew Knepley 19 Aug '05

19 Aug '05

Simon Burton <simon(a)arrowtheory.com> writes: > On Thu, 18 Aug 2005 19:57:02 -0500 > Matthew Knepley <knepley(a)mcs.anl.gov> wrote: > >> >> > (d) m=exp(m) (pointwise exp) >> >> We could add pointwise operations just like the VecPointwise*(). >> >> Matt > > Yes, we need MatPointwiseMult aswell. > > Should I make a start on this ? I am thinking more about what Barry said. The VecPointwise*() operations can be given a solid mathematical interpretation in terms of spinor operations. However, I do not see anything like that for the Mat stuff yet. We need to understand the mathematicas better. Matt -- "Failure has a thousand explanations. Success doesn't need one" -- Sir Alec Guiness

1 0

Re: ufuncs, iterators
by Matthew Knepley 19 Aug '05

19 Aug '05

S V N Vishwanathan <vishy(a)mail.rsise.anu.edu.au> writes: > Hi! > >>> BTW: possibly related note, are you using dense matrices sometimes >>> to represent "just" 2-arrays; that is, not as representations of linear >>> operators. If so, I do not think this is the correct approach! Conceptually >>> PETSc Mat's are linear operators I think it would be a big mistake to >>> overload them as 2-arrays also. The "correct" approach is to use the >>> DACreate2d() construct for handling 2-d arrays; With the DA the values >>> are stored into Vecs but there is additional information about >>> the two array structure; it can be decomposed nicely in parallel and >>> one can set/access values with the usual two i,j indices. Of course if >>> they are being used as operators ignore this. > > Simon> OK, I will think about this. But what would be the strategy for > Simon> doing all these linear algebra operations ? Swap to/from Mat/DA > Simon> objects ? > >>> In Matlab matrix == 2-d array == matrix (in my opinion a terrible >>> design decision) in PETSc matrix (dense) != 2-d array , the are >>> completely different beasts mathematically > > I am not sure I understand the fine difference. As far as we are > concerned, all the operations which we are doing (point wise addition, > addition, multiplication etc.) are on the linear operator. But it might > be that my thought process is conditioned by years of Matlab/Octave > use. Can you maybe make this more explicit i.e. which situations would > you use a DA array and when would you use a Mat object? Linear operators are members of the space of linear maps between two vector spaces. They, for instance, have a characteristic behavior under coordinate transformations. 2-d arrays are just lists of values. They need not have any transformation properties. In fact, linear operators also have spectral characteristics, and analytic behaviors. Using a PETSc Mat just to store a collection of values is wrong. This goes back to the modeling for your equations. Matt -- "Failure has a thousand explanations. Success doesn't need one" -- Sir Alec Guiness

1 0

Re: ufuncs, iterators
by Barry Smith 18 Aug '05

18 Aug '05

In terms of exp( -1/2\sigma^{2} ||x_{i} - x_{j}||_{2}^{2}) what are they? Thanks Barry On Fri, 19 Aug 2005, Simon Burton wrote: > On Thu, 18 Aug 2005 22:49:29 -0500 (CDT) > Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > > > > What is x1, x2 and ip? > > > > Barry > > x1 and x2 are 2-arrays; their rows are the 'sample' vectors. > ip is the matrix of all inner products from x1 and x2. > > Simon. > > >

1 0

Re: ufuncs, iterators
by Simon Burton 18 Aug '05

18 Aug '05

On Thu, 18 Aug 2005 22:49:29 -0500 (CDT) Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > What is x1, x2 and ip? > > Barry x1 and x2 are 2-arrays; their rows are the 'sample' vectors. ip is the matrix of all inner products from x1 and x2. Simon. -- Simon Burton, B.Sc. Licensed PO Box 8066 ANU Canberra 2601 Australia Ph. 61 02 6249 6940 http://arrowtheory.com

1 0

Re: ufuncs, iterators
by Barry Smith 18 Aug '05

18 Aug '05

What is x1, x2 and ip? Barry On Fri, 19 Aug 2005, S V N Vishwanathan wrote: > Hi! > > >> I am not sure I understand the fine difference. As far as we are > >> concerned, all the operations which we are doing (point wise addition, > >> addition, multiplication etc.) are on the linear operator. > > Barry> Certainly pointwise addition is the same as adding two operators > Barry> together (PETSc has this with MatAXPY), but what about, for > Barry> example, pointwise multiply? You could well be correct, but I'd > Barry> like to see the list and what they correspond to in terms of > Barry> linear operators, for example row and column sums what do they > Barry> represent (true they are the multiplication of the matrix (or its > Barry> transpose) by the vector of all ones, but is that a useful > Barry> meaning. > > Well I guess the problem is that we are interested in computing matrices > of the form K_{ij} where the i, j th entry is > > exp( -1/2\sigma^{2} ||x_{i} - x_{j}||_{2}^{2}) > > Here X is a huge matrix (possibly sparse) and x_{i} are the rows of > X. Of course, in a purist sense this is not linear algebra being applied > to a linear operator but it makes sense to use pointwise operations to > compute it. No? > > Sample code: > > n_sq1 = numarray.add.reduce( x1*x1, 1 ) # vector of norm squareds (slightly faster than x1**2) > n_sq2 = numarray.add.reduce( x2*x2, 1 ) > n_sq = (-0.5/self.sigma2) * numarray.add.outer( n_sq1, n_sq2 ) # matrix of norm squared sums > res = ip # result (alias) > numarray.add( ip, n_sq, res ) > numarray.exp( res, res ) > > vishy > >

1 0

Re: ufuncs, iterators
by S V N Vishwanathan 18 Aug '05

18 Aug '05

Hi! >> I am not sure I understand the fine difference. As far as we are >> concerned, all the operations which we are doing (point wise addition, >> addition, multiplication etc.) are on the linear operator. Barry> Certainly pointwise addition is the same as adding two operators Barry> together (PETSc has this with MatAXPY), but what about, for Barry> example, pointwise multiply? You could well be correct, but I'd Barry> like to see the list and what they correspond to in terms of Barry> linear operators, for example row and column sums what do they Barry> represent (true they are the multiplication of the matrix (or its Barry> transpose) by the vector of all ones, but is that a useful Barry> meaning. Well I guess the problem is that we are interested in computing matrices of the form K_{ij} where the i, j th entry is exp( -1/2\sigma^{2} ||x_{i} - x_{j}||_{2}^{2}) Here X is a huge matrix (possibly sparse) and x_{i} are the rows of X. Of course, in a purist sense this is not linear algebra being applied to a linear operator but it makes sense to use pointwise operations to compute it. No? Sample code: n_sq1 = numarray.add.reduce( x1*x1, 1 ) # vector of norm squareds (slightly faster than x1**2) n_sq2 = numarray.add.reduce( x2*x2, 1 ) n_sq = (-0.5/self.sigma2) * numarray.add.outer( n_sq1, n_sq2 ) # matrix of norm squared sums res = ip # result (alias) numarray.add( ip, n_sq, res ) numarray.exp( res, res ) vishy

1 0

Re: ufuncs, iterators
by Barry Smith 18 Aug '05

18 Aug '05

On Fri, 19 Aug 2005, S V N Vishwanathan wrote: > Hi! > > >> BTW: possibly related note, are you using dense matrices sometimes > >> to represent "just" 2-arrays; that is, not as representations of linear > >> operators. If so, I do not think this is the correct approach! Conceptually > >> PETSc Mat's are linear operators I think it would be a big mistake to > >> overload them as 2-arrays also. The "correct" approach is to use the > >> DACreate2d() construct for handling 2-d arrays; With the DA the values > >> are stored into Vecs but there is additional information about > >> the two array structure; it can be decomposed nicely in parallel and > >> one can set/access values with the usual two i,j indices. Of course if > >> they are being used as operators ignore this. > > Simon> OK, I will think about this. But what would be the strategy for > Simon> doing all these linear algebra operations ? Swap to/from Mat/DA > Simon> objects ? > > >> In Matlab matrix == 2-d array == matrix (in my opinion a terrible > >> design decision) in PETSc matrix (dense) != 2-d array , the are > >> completely different beasts mathematically > > I am not sure I understand the fine difference. As far as we are > concerned, all the operations which we are doing (point wise addition, > addition, multiplication etc.) are on the linear operator. Certainly pointwise addition is the same as adding two operators together (PETSc has this with MatAXPY), but what about, for example, pointwise multiply? You could well be correct, but I'd like to see the list and what they correspond to in terms of linear operators, for example row and column sums what do they represent (true they are the multiplication of the matrix (or its transpose) by the vector of all ones, but is that a useful meaning. I just want to make sure we are using the right paradigm for these operations, before we move to far along in where we implement them. Barry > But it might > be that my thought process is conditioned by years of Matlab/Octave > use. Can you maybe make this more explicit i.e. which situations would > you use a DA array and when would you use a Mat object? > > vishy > >

1 0

Re: ufuncs, iterators
by S V N Vishwanathan 18 Aug '05

18 Aug '05

Hi! >> BTW: possibly related note, are you using dense matrices sometimes >> to represent "just" 2-arrays; that is, not as representations of linear >> operators. If so, I do not think this is the correct approach! Conceptually >> PETSc Mat's are linear operators I think it would be a big mistake to >> overload them as 2-arrays also. The "correct" approach is to use the >> DACreate2d() construct for handling 2-d arrays; With the DA the values >> are stored into Vecs but there is additional information about >> the two array structure; it can be decomposed nicely in parallel and >> one can set/access values with the usual two i,j indices. Of course if >> they are being used as operators ignore this. Simon> OK, I will think about this. But what would be the strategy for Simon> doing all these linear algebra operations ? Swap to/from Mat/DA Simon> objects ? >> In Matlab matrix == 2-d array == matrix (in my opinion a terrible >> design decision) in PETSc matrix (dense) != 2-d array , the are >> completely different beasts mathematically I am not sure I understand the fine difference. As far as we are concerned, all the operations which we are doing (point wise addition, addition, multiplication etc.) are on the linear operator. But it might be that my thought process is conditioned by years of Matlab/Octave use. Can you maybe make this more explicit i.e. which situations would you use a DA array and when would you use a Mat object? vishy

1 0

Re: ufuncs, iterators
by Simon Burton 18 Aug '05

18 Aug '05

On Thu, 18 Aug 2005 20:50:12 -0500 (CDT) Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > Simon, > > There is not because what does MatSet() mean for a sparse matrix? > Make it dense, just set all the current nonzeros, error? You tell us, > I could live with either of the later two but don't really like the > first. Yes, it would be an error, that sounds the best to me. > > Barry > > BTW: possibly related note, are you using dense matrices sometimes > to represent "just" 2-arrays; that is, not as representations of linear > operators. If so, I do not think this is the correct approach! Conceptually > PETSc Mat's are linear operators I think it would be a big mistake to > overload them as 2-arrays also. The "correct" approach is to use the > DACreate2d() construct for handling 2-d arrays; With the DA the values > are stored into Vecs but there is additional information about > the two array structure; it can be decomposed nicely in parallel and > one can set/access values with the usual two i,j indices. Of course if > they are being used as operators ignore this. OK, I will think about this. But what would be the strategy for doing all these linear algebra operations ? Swap to/from Mat/DA objects ? Simon. > > In Matlab matrix == 2-d array == matrix (in my opinion a terrible design > decision) > in PETSc matrix (dense) != 2-d array , the are completely different > beasts mathematically > > > > On Fri, 19 Aug 2005, Simon Burton wrote: > > > > > Is there a way to simulate a MatSet (like VecSet) ? > > > > Simon. > > > > -- Simon Burton, B.Sc. Licensed PO Box 8066 ANU Canberra 2601 Australia Ph. 61 02 6249 6940 http://arrowtheory.com

1 0