Hi all,
-- I've been looking at DhbVector implementation and couple of questions came up:
How I would like to see it:
I am not sure if this is achievable and what kind of design should be chosen to implement this with the minimum code (once again, in Common Lisp I would probably resort to macros as the code is essentially the same and only type hints are changing). On the top of this, to efficiently implement some extra linear algebra methods more BLAS-like operations on vectors are needed, i.e. operations that don't allocate the memory but modify existing vector. This can be used in, for example, implementation of Conjugate Gradients (CG) linear solver for positively defined symmetric matrices, BiCGStab and GEMRES for general matrices; these methods are efficient for sparse systems and can solve linear equations where matrix part is expressed as a function. I'm keen to add these operations: I can start adding them to DhbVector and if we decide to change the whole vector design, they can be moved to an appropriate class. Cheers, Alexey PS: scmutils, the extension of MIT-Scheme, to study Classical Mechanics developed by G. Sussman, A. Radul and others from MIT uses an interesting idea for vectors: there are two kinds of vectors, columns (vectors, with indices in superscripts) and rows (co-vectors, with indices in subscripts); so, scalar multiplication would be achieved by multiplying a row by a column; if a column is multiplied by a row, this will give outer product. Furthermore, element-wise operations are only allowed for the same kind of vectors. You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
On Monday, April 4, 2016 at 4:10:45 PM UTC+2, Alexey Cherkaev wrote:
You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
In reply to this post by Alexey Cherkaev
Hi Alexey,
-- calculations with Floatarrays are of course faster but the precision is much worse than than normal arrays, imo it would be not so good to loose so much precision. werner On Monday, April 4, 2016 at 4:10:45 PM UTC+2, Alexey Cherkaev wrote:
You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
example:
-- a:=Array with:Float pi. f:=a asFloatArray . a squared. "#(9.869604401089358)" f squared. "a FloatArray(9.86960506439209)" You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
In reply to this post by Alexey Cherkaev
The class DHBVector has been introduced to implement common linear algebra operations: product with scalar, addition, dot product and matrix products.Using float is not good enough for precision as has laready been addressed. As you pointed out, there was no need for additional instance variables, only additonal methods and operators were added.Those methods are needed in all linear algebra algorithms. I would also add that the methods and operators have been designed to improve readability. For further explanation, I refer you to my book: OO implementation of numerical methods. Cheers, Didier On Mon, Apr 4, 2016 at 4:10 PM, Alexey Cherkaev <[hidden email]> wrote:
You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
Hi Werner and Didier,
-- Thanks for the replies! I didn't expect the precision deteriorate in FloatArray... Do you know what is the reason? So, my understanding of DhbVector is that it represents a) a mathematical vector (conceptually); b) an array of some kind of numbers, including Float, DualNumber, etc.). It is a different entity from Array, as Array is a more generic collection of things (ordered and indexed). Am I correct? (my question if the separate class was needed was based on the idea of using FloatArray, before realising it was losing the precision) As for extra operations I find myself using (while implementing CG method): y <- a * y (1) y <- y + a * x (2) y <- a * y + x (3) y <- y + x (4) with x, y - vectors and a - scalar. And (1) is already in: scaleBy method. On Monday, April 4, 2016 at 9:26:35 PM UTC+2, Didier Besset wrote:
You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
2016-04-04 22:46 GMT+02:00 Alexey Cherkaev <[hidden email]>:
FloatArray is single precision (32 bits) floats While Squeak/Pharo Float is double precision (64 bits)
This sounds too much like BLAS primitives...
You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
In reply to this post by Alexey Cherkaev
Hi,
-- re #scaleBy: long time ago i made #* commutatively usable with any mixture of scalars, dhbvectors and dhbmatrices and replaced all uses of #scaleBy: in dhb by #*. i thought using #* is conceptually simpler and one could perhaps one day deprecate #scaleBy: or so. werner You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
oops, of course i mean the mixture of scalars with dhbvectors and dhbmatrices are commutatively mixable.
-- On Tuesday, April 5, 2016 at 10:57:55 AM UTC+2, werner kassens wrote:
You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
In reply to this post by werner kassens-2
Hi Nicolas, Werner!
-- I must admit, Pharo's non-distinction on the naming level between double-float for Float numbers and single-float FloatArray took me by surprise. Werner: #scaleBy: does not produce a new vector as a result, but modifies the receiver. Whereas, I assume, it is expected that 'aNumber * aVector' doesn't modify aVector and answers a newly allocated vector. #scaleBy: in that sense is a lower level operation. Which brings me to the next point: Nicolas: yes, these operations are primitive. The reason for them is that they are used repeatedly in iterative solvers (CG, BiCGStab) which can be used in non-linear solvers (e.g. Newton's method), which can be used in ODE sovlers.... You get the picture. But, I was thinking about those primitives: they are not really useful to the user of the code. The user would want to operate using #+, #* and wouldn't really care about producing (and trashing) some intermediate results. So, here is my dilemma: I would want to make CG, BiCGStab and GEMRES (this one is a bit trickier) efficient, so that they can be used inside non-linear solvers. While I've been playing with CG and BiCGStab in Common Lisp, I found that the solution timing is dominated by matrix to vector multiplication (or equivalently, linear operator to vector application), even for sparse systems. So, actual vector arithmetic is fast. Yet, I'm afraid that producing lots of intermediate short-lived results will strain garbage collector, once these methods become a part of non-linear solver. Should I just implement these methods using existing vector arithmetic and worry about optimisation later? I.e. make it work, make it right, make it fast approach. (I do feel the answer to this question should be 'yes'). Secondly, in terms of the future of PolyMath, should DhbVector include BLAS-like primitive? It might feel like over-engineering at this point though. Nice thing about these primitives, however, is that higher-level operations can be easily implemented using these primitives. So, it introduces nice level of abstraction. (and, let's say in the future, these primitives can be implemented with parallelisation on GPU) Thanks for the help! Alexey On Tuesday, April 5, 2016 at 10:57:55 AM UTC+2, werner kassens wrote:
You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
On Tuesday, April 5, 2016 at 11:52:47 AM UTC+2, Alexey Cherkaev wrote:
--
Hi Alexey, yes, that is why i did not silently deprecate #scaleBy: <stupid grin>. in dhb Didier often uses in place calcs instead of producing intermediary objects, that then gets thrown away, but he often put them in a bigger method where several of these object updates are done in one iteration. wouldnt this be a possibility instead of a method that just updates y < - y + a * x etc? werner You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
Hi Werner,
-- Not sure if I fully understand what you've meant. Can you please elaborate a bit? Also, if you look at CG algorithm https://en.wikipedia.org/wiki/Conjugate_gradient_method#The_resulting_algorithm, apart from updating the solution and residual, there are no multiple vector updates, that wouldn't be punctuated by something like scalar product. But, perhaps, x = x + alpha p, r = r - alpha (A p) and beta = (r-new, r-new) / (r,-old, r-old) can be done in one iteration ((A p) would be precomputed t that point). Cheers, Alexey On Tuesday, April 5, 2016 at 12:19:23 PM UTC+2, werner kassens wrote:
You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
In reply to this post by Alexey Cherkaev
2016-04-05 11:52 GMT+02:00 Alexey Cherkaev <[hidden email]>:
BLAS certainly makes sense if you want that basic linear algebra operations on numerical data to better scale (in terms of operations/second especially for large matrices, because it's not really sure it applies to 3x3 dimensions) If we want to compete with, say a basic Matlab interpreter or Numpy or R or etc..., then it's mandatory from my POV. If we are using a BLAS backend, then an idea is to create proxies and don't perform any operation untill it better matches a native BLAS operation. This is implemented in C++ for example in BOOST ublas http://www.boost.org/doc/libs/1_60_0/libs/numeric/ublas/doc/ May I remind that I have developped a rather rough and dumb interface to BLACK/LAPACK for Smalltalk: Smallapack see https://github.com/nicolas-cellier-aka-nice/smallapack and https://github.com/nicolas-cellier-aka-nice/smallapack/wiki Its status on latest Pharo is unknown, the package depends on old Compiler, but if there's interest into it I can inquire. Smallapack is not integrated in PolyMath because maybe it's too huge, and because interaction with dhb matrices/vector would have to be thought out first. I see it as a concurrent implementation, so one could choose either dhb or blas/lapack based matrix/vector implementation depending on the needs, and we could maintain a kind of minimal common API between the two. Nicolas
You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
As a disclaimer: I am definitely not in a position to make this call, I guess Stephan, Serge and Didier would need to decide. I can comment as a user: my potential use of PolyMath would include solution of ODE of size ~100, with minimization problem solution on the top of it. So, on one hand I don't need huge scalability. On the other hand, this problem requires an efficient foundation and BLAS/LAPACK may provide it. PS. What is the current state of FFI in Pharo? I saw there was a talk on this topic at PharoDays. You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
In reply to this post by Nicolas Cellier
On Tue, Apr 5, 2016 at 5:44 PM, Nicolas Cellier
<[hidden email]> wrote: > > > 2016-04-05 11:52 GMT+02:00 Alexey Cherkaev <[hidden email]>: >> >> Hi Nicolas, Werner! >> >> I must admit, Pharo's non-distinction on the naming level between >> double-float for Float numbers and single-float FloatArray took me by >> surprise. >> >> Werner: #scaleBy: does not produce a new vector as a result, but modifies >> the receiver. Whereas, I assume, it is expected that 'aNumber * aVector' >> doesn't modify aVector and answers a newly allocated vector. #scaleBy: in >> that sense is a lower level operation. Which brings me to the next point: >> >> Nicolas: yes, these operations are primitive. The reason for them is that >> they are used repeatedly in iterative solvers (CG, BiCGStab) which can be >> used in non-linear solvers (e.g. Newton's method), which can be used in ODE >> sovlers.... You get the picture. But, I was thinking about those primitives: >> they are not really useful to the user of the code. The user would want to >> operate using #+, #* and wouldn't really care about producing (and trashing) >> some intermediate results. >> >> So, here is my dilemma: >> >> I would want to make CG, BiCGStab and GEMRES (this one is a bit trickier) >> efficient, so that they can be used inside non-linear solvers. While I've >> been playing with CG and BiCGStab in Common Lisp, I found that the solution >> timing is dominated by matrix to vector multiplication (or equivalently, >> linear operator to vector application), even for sparse systems. So, actual >> vector arithmetic is fast. Yet, I'm afraid that producing lots of >> intermediate short-lived results will strain garbage collector, once these >> methods become a part of non-linear solver. >> >> Should I just implement these methods using existing vector arithmetic and >> worry about optimisation later? I.e. make it work, make it right, make it >> fast approach. (I do feel the answer to this question should be 'yes'). >> >> Secondly, in terms of the future of PolyMath, should DhbVector include >> BLAS-like primitive? It might feel like over-engineering at this point >> though. Nice thing about these primitives, however, is that higher-level >> operations can be easily implemented using these primitives. So, it >> introduces nice level of abstraction. (and, let's say in the future, these >> primitives can be implemented with parallelisation on GPU) >> >> Thanks for the help! >> Alexey >> > > BLAS certainly makes sense if you want that basic linear algebra operations > on numerical data to better scale > (in terms of operations/second especially for large matrices, because it's > not really sure it applies to 3x3 dimensions) > If we want to compete with, say a basic Matlab interpreter or Numpy or R or > etc..., then it's mandatory from my POV. > > If we are using a BLAS backend, then an idea is to create proxies and don't > perform any operation untill it better matches a native BLAS operation. > This is implemented in C++ for example in BOOST ublas > http://www.boost.org/doc/libs/1_60_0/libs/numeric/ublas/doc/ > > May I remind that I have developped a rather rough and dumb interface to > BLACK/LAPACK for Smalltalk: Smallapack > see https://github.com/nicolas-cellier-aka-nice/smallapack and > https://github.com/nicolas-cellier-aka-nice/smallapack/wiki > Its status on latest Pharo is unknown, the package depends on old Compiler, > but if there's interest into it I can inquire. > > Smallapack is not integrated in PolyMath because maybe it's too huge, and > because interaction with dhb matrices/vector would have to be thought out > first. > I see it as a concurrent implementation, so one could choose either dhb or > blas/lapack based matrix/vector implementation depending on the needs, and > we could maintain a kind of minimal common API between the two. Would be great to integrate your work on Smallapack later in PolyMath ! At the moment, we have to finish to clean the code, the dependancies between the packages and uniform prefix. Yes you are right, I guess we could try to define a common API interface. Do you have time to work on this ? -- Serge Stinckwich UCBN & UMI UMMISCO 209 (IRD/UPMC) Every DSL ends up being Smalltalk http://www.doesnotunderstand.org/ -- You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
In reply to this post by Alexey Cherkaev
On Tue, Apr 5, 2016 at 6:36 PM, Alexey Cherkaev
<[hidden email]> wrote: >> BLAS certainly makes sense if you want that basic linear algebra >> operations on numerical data to better scale >> (in terms of operations/second especially for large matrices, because it's >> not really sure it applies to 3x3 dimensions) >> If we want to compete with, say a basic Matlab interpreter or Numpy or R >> or etc..., then it's mandatory from my POV. >> >> If we are using a BLAS backend, then an idea is to create proxies and >> don't perform any operation untill it better matches a native BLAS >> operation. >> This is implemented in C++ for example in BOOST ublas >> http://www.boost.org/doc/libs/1_60_0/libs/numeric/ublas/doc/ >> >> May I remind that I have developped a rather rough and dumb interface to >> BLACK/LAPACK for Smalltalk: Smallapack >> see https://github.com/nicolas-cellier-aka-nice/smallapack and >> https://github.com/nicolas-cellier-aka-nice/smallapack/wiki >> Its status on latest Pharo is unknown, the package depends on old >> Compiler, but if there's interest into it I can inquire. >> >> Smallapack is not integrated in PolyMath because maybe it's too huge, and >> because interaction with dhb matrices/vector would have to be thought out >> first. >> I see it as a concurrent implementation, so one could choose either dhb or >> blas/lapack based matrix/vector implementation depending on the needs, and >> we could maintain a kind of minimal common API between the two. >> >> Nicolas > > > As a disclaimer: I am definitely not in a position to make this call, I > guess Stephan, Serge and Didier would need to decide. I can comment as a > user: my potential use of PolyMath would include solution of ODE of size > ~100, with minimization problem solution on the top of it. So, on one hand I > don't need huge scalability. On the other hand, this problem requires an > efficient foundation and BLAS/LAPACK may provide it. This will be really great if we can have problems with size that run on PolyMath. We are user-driven here, please propose code and scenarios to do in that direction :-) -- Serge Stinckwich UCBN & UMI UMMISCO 209 (IRD/UPMC) Every DSL ends up being Smalltalk http://www.doesnotunderstand.org/ -- You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
In reply to this post by Alexey Cherkaev
Hi Alexey,
-- i have to admit i could not exactly identify your mentioned updates (1 to 4) in that algo, but i assume we are talking about x(i), p(i), r(i) ? they have all the same size, correct? in that case one could eg try for optimizing purposes to do a 1to:size do[:position| x(i) at:position put:whatever.p(i)at:position bla.r(i)bla] or something like that. (and one would probably do the a(i)*p(i) beforehand and perhaps some other minor points, i had just a short glance at this algo. the matrix multiplication could also be done stepwise in this thing, or? ). is that what you have in mind, or do i misunderstand you? werner You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
In reply to this post by Nicolas Cellier
On Tuesday, April 5, 2016 at 5:44:07 PM UTC+2, Nicolas Cellier wrote:
--
Hi Nicolas, i always feared that smallapack is to complicated to use for me, hence i never tried it. a minimal common API would probably help me overcoming that fear. (and i would need some simple installation instructions) werner You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
In reply to this post by werner kassens-2
i just noticed one cant do it in one go, beta has to be calced after the first go and in a second one p(i).
-- On Tuesday, April 5, 2016 at 7:38:33 PM UTC+2, werner kassens wrote:
You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
In reply to this post by Nicolas Cellier
Hello nicolas
if you want we could add/include Smallapack into polymath. In the long term I would like to split the huge configurations in many smaller one and the user should be able to pick what he wants Now I'm on holidays trying to get far from a keyboard. Stef Le 5/4/16 17:44, Nicolas Cellier a
écrit :
-- You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
Free forum by Nabble - Scala forum | Edit this page |