Smalltalk › Frameworks & Tools › SciSmalltalk

DhbVector: implementation and extensions

_‹ Previous Topic Next Topic _›

Classic

List

Threaded

27 messages Options

Alexey Cherkaev

DhbVector: implementation and extensions

Hi all,

I've been looking at DhbVector implementation and couple of questions came up:

DhbVector is a sublcass of Array. Why not FloatArray? I see PolyMath has a package for automatic differentiation. If AD is used, the user will need the array of DualNumber's or Tape's (not implemented yet) In this case FloatArray might get on the way.
Even more, why we even need a separate class for a vector? DhbVector doesn't have any instance or class variable. Vector operations can be implemented as extensions of Array/FloatArray.

How I would like to see it:

Have FloatArray implementing all the necessary operations for Float's only. The idea is that this would be an optimised implementation. In Common Lisp I would implement operations as multi-methods with type hints inside specialisations for floating-point arrays. I don't know how to achieve this in Pharo and what kind of optimisation it can offer.
Have some kind of NumberArray, whose members are known to be some kind of numbers. So, it will be able to deal with DualNumber's, for example. Thus, it's less optimised, but more universal.

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.

werner kassens-2

Re: DhbVector: implementation and extensions

On Monday, April 4, 2016 at 4:10:45 PM UTC+2, Alexey Cherkaev wrote:

Hi all,

I've been looking at DhbVector implementation and couple of questions came up:
DhbVector is a sublcass of Array. Why not FloatArray? I see PolyMath has a package for automatic differentiation. If AD is used, the user will need the array of DualNumber's or Tape's (not implemented yet) In this case FloatArray might get on the way.
Even more, why we even need a separate class for a vector? DhbVector doesn't have any instance or class variable. Vector operations can be implemented as extensions of Array/FloatArray.
How I would like to see it:
Have FloatArray implementing all the necessary operations for Float's only. The idea is that this would be an optimised implementation. In Common Lisp I would implement operations as multi-methods with type hints inside specialisations for floating-point arrays. I don't know how to achieve this in Pharo and what kind of optimisation it can offer.
Have some kind of NumberArray, whose members are known to be some kind of numbers. So, it will be able to deal with DualNumber's, for example. Thus, it's less optimised, but more universal.
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.

werner kassens-2

Re: DhbVector: implementation and extensions

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:

Hi all,

I've been looking at DhbVector implementation and couple of questions came up:
DhbVector is a sublcass of Array. Why not FloatArray? I see PolyMath has a package for automatic differentiation. If AD is used, the user will need the array of DualNumber's or Tape's (not implemented yet) In this case FloatArray might get on the way.
Even more, why we even need a separate class for a vector? DhbVector doesn't have any instance or class variable. Vector operations can be implemented as extensions of Array/FloatArray.
How I would like to see it:
Have FloatArray implementing all the necessary operations for Float's only. The idea is that this would be an optimised implementation. In Common Lisp I would implement operations as multi-methods with type hints inside specialisations for floating-point arrays. I don't know how to achieve this in Pharo and what kind of optimisation it can offer.
Have some kind of NumberArray, whose members are known to be some kind of numbers. So, it will be able to deal with DualNumber's, for example. Thus, it's less optimised, but more universal.
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.

werner kassens-2

Re: DhbVector: implementation and extensions

example:
a:=Array with:Float pi.
f:=a asFloatArray .
a squared. "#(9.869604401089358)"
f squared. "a FloatArray(9.86960506439209)"

Didier Besset

Re: DhbVector: implementation and extensions

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:

Hi all,

I've been looking at DhbVector implementation and couple of questions came up:
DhbVector is a sublcass of Array. Why not FloatArray? I see PolyMath has a package for automatic differentiation. If AD is used, the user will need the array of DualNumber's or Tape's (not implemented yet) In this case FloatArray might get on the way.
Even more, why we even need a separate class for a vector? DhbVector doesn't have any instance or class variable. Vector operations can be implemented as extensions of Array/FloatArray.
How I would like to see it:
Have FloatArray implementing all the necessary operations for Float's only. The idea is that this would be an optimised implementation. In Common Lisp I would implement operations as multi-methods with type hints inside specialisations for floating-point arrays. I don't know how to achieve this in Pharo and what kind of optimisation it can offer.
Have some kind of NumberArray, whose members are known to be some kind of numbers. So, it will be able to deal with DualNumber's, for example. Thus, it's less optimised, but more universal.
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.

Alexey Cherkaev

Re: DhbVector: implementation and extensions

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:

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 <<a href="javascript:" target="_blank" gdf-obfuscated-mailto="oK9AtnMiCAAJ" rel="nofollow" onmousedown="this.href='javascript:';return true;" onclick="this.href='javascript:';return true;">alexey....@...> wrote:
Hi all,

I've been looking at DhbVector implementation and couple of questions came up:
DhbVector is a sublcass of Array. Why not FloatArray? I see PolyMath has a package for automatic differentiation. If AD is used, the user will need the array of DualNumber's or Tape's (not implemented yet) In this case FloatArray might get on the way.
Even more, why we even need a separate class for a vector? DhbVector doesn't have any instance or class variable. Vector operations can be implemented as extensions of Array/FloatArray.
How I would like to see it:
Have FloatArray implementing all the necessary operations for Float's only. The idea is that this would be an optimised implementation. In Common Lisp I would implement operations as multi-methods with type hints inside specialisations for floating-point arrays. I don't know how to achieve this in Pharo and what kind of optimisation it can offer.
Have some kind of NumberArray, whose members are known to be some kind of numbers. So, it will be able to deal with DualNumber's, for example. Thus, it's less optimised, but more universal.
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 <a href="javascript:" target="_blank" gdf-obfuscated-mailto="oK9AtnMiCAAJ" rel="nofollow" onmousedown="this.href='javascript:';return true;" onclick="this.href='javascript:';return true;">scismalltalk...@googlegroups.com.
For more options, visit <a href="https://groups.google.com/d/optout" target="_blank" rel="nofollow" onmousedown="this.href='https://groups.google.com/d/optout';return true;" onclick="this.href='https://groups.google.com/d/optout';return true;">https://groups.google.com/d/optout.

Nicolas Cellier

Re: DhbVector: implementation and extensions

2016-04-04 22:46 GMT+02:00 Alexey Cherkaev <[hidden email]>:

Hi Werner and Didier,

Thanks for the replies! I didn't expect the precision deteriorate in FloatArray... Do you know what is the reason?

FloatArray is single precision (32 bits) floats

While Squeak/Pharo Float is double precision (64 bits)

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.

This sounds too much like BLAS primitives...

On Monday, April 4, 2016 at 9:26:35 PM UTC+2, Didier Besset wrote:
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:
Hi all,

I've been looking at DhbVector implementation and couple of questions came up:
DhbVector is a sublcass of Array. Why not FloatArray? I see PolyMath has a package for automatic differentiation. If AD is used, the user will need the array of DualNumber's or Tape's (not implemented yet) In this case FloatArray might get on the way.
Even more, why we even need a separate class for a vector? DhbVector doesn't have any instance or class variable. Vector operations can be implemented as extensions of Array/FloatArray.
How I would like to see it:
Have FloatArray implementing all the necessary operations for Float's only. The idea is that this would be an optimised implementation. In Common Lisp I would implement operations as multi-methods with type hints inside specialisations for floating-point arrays. I don't know how to achieve this in Pharo and what kind of optimisation it can offer.
Have some kind of NumberArray, whose members are known to be some kind of numbers. So, it will be able to deal with DualNumber's, for example. Thus, it's less optimised, but more universal.
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.

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

werner kassens-2

Re: DhbVector: implementation and extensions

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

werner kassens-2

Re: DhbVector: implementation and extensions

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:

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

Alexey Cherkaev

Re: DhbVector: implementation and extensions

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:

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

werner kassens-2

Re: DhbVector: implementation and extensions

On Tuesday, April 5, 2016 at 11:52:47 AM UTC+2, Alexey Cherkaev wrote:

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

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

Alexey Cherkaev

Re: DhbVector: implementation and extensions

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:

On Tuesday, April 5, 2016 at 11:52:47 AM UTC+2, Alexey Cherkaev wrote:
#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:

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

Nicolas Cellier

Re: DhbVector: implementation and extensions

In reply to this post by Alexey Cherkaev

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.

Nicolas

On Tuesday, April 5, 2016 at 10:57:55 AM UTC+2, werner kassens wrote:
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.

Alexey Cherkaev

Re: DhbVector: implementation and extensions

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 <a href="http://www.boost.org/doc/libs/1_60_0/libs/numeric/ublas/doc/" target="_blank" rel="nofollow" onmousedown="this.href='http://www.google.com/url?q\75http%3A%2F%2Fwww.boost.org%2Fdoc%2Flibs%2F1_60_0%2Flibs%2Fnumeric%2Fublas%2Fdoc%2F\46sa\75D\46sntz\0751\46usg\75AFQjCNFUfHEzMVzfo_uBjJh340Des8ulgg';return true;" onclick="this.href='http://www.google.com/url?q\75http%3A%2F%2Fwww.boost.org%2Fdoc%2Flibs%2F1_60_0%2Flibs%2Fnumeric%2Fublas%2Fdoc%2F\46sa\75D\46sntz\0751\46usg\75AFQjCNFUfHEzMVzfo_uBjJh340Des8ulgg';return true;">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 <a href="https://github.com/nicolas-cellier-aka-nice/smallapack" target="_blank" rel="nofollow" onmousedown="this.href='https://www.google.com/url?q\75https%3A%2F%2Fgithub.com%2Fnicolas-cellier-aka-nice%2Fsmallapack\46sa\75D\46sntz\0751\46usg\75AFQjCNFmQmaktlrBIAmVuylhO1luF04f-Q';return true;" onclick="this.href='https://www.google.com/url?q\75https%3A%2F%2Fgithub.com%2Fnicolas-cellier-aka-nice%2Fsmallapack\46sa\75D\46sntz\0751\46usg\75AFQjCNFmQmaktlrBIAmVuylhO1luF04f-Q';return true;">https://github.com/nicolas-cellier-aka-nice/smallapack and <a href="https://github.com/nicolas-cellier-aka-nice/smallapack/wiki" target="_blank" rel="nofollow" onmousedown="this.href='https://www.google.com/url?q\75https%3A%2F%2Fgithub.com%2Fnicolas-cellier-aka-nice%2Fsmallapack%2Fwiki\46sa\75D\46sntz\0751\46usg\75AFQjCNGOiHE6FFZWcuTLyrUAERTSdt7D-Q';return true;" onclick="this.href='https://www.google.com/url?q\75https%3A%2F%2Fgithub.com%2Fnicolas-cellier-aka-nice%2Fsmallapack%2Fwiki\46sa\75D\46sntz\0751\46usg\75AFQjCNGOiHE6FFZWcuTLyrUAERTSdt7D-Q';return true;">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.

PS. What is the current state of FFI in Pharo? I saw there was a talk on this topic at PharoDays.

SergeStinckwich

Re: DhbVector: implementation and extensions

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.

SergeStinckwich

Re: DhbVector: implementation and extensions

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.

werner kassens-2

Re: DhbVector: implementation and extensions

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

werner kassens-2

Re: DhbVector: implementation and extensions

In reply to this post by Nicolas Cellier

On Tuesday, April 5, 2016 at 5:44:07 PM UTC+2, Nicolas Cellier wrote:

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.

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

werner kassens-2

Re: DhbVector: implementation and extensions

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:

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

stepharo

Re: DhbVector: implementation and extensions

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 :

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.

Nicolas

On Tuesday, April 5, 2016 at 10:57:55 AM UTC+2, werner kassens wrote:

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.

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