From: Jeronimo Pellegrini
Subject: Sparse matrices in Common Lisp? (with 3 dimentions)
Date:
Message-ID: <f3qi3f$vp1$1@aioe.org>
Hello.
I was wondering what would be an efficient way to implement
three-dimensional sparse matrices in Common Lisp (from what I see in the
HyperSpec, there is no native support for them).
I also need them to have fast access...
I thought I would implement them as hash-tables, with wrappers over
these two:
(setf (gethash (list i j k) my-matrix) value)
(gethash (list i j k) my-matrix)
And then have some "commit" function that will map the non-zeroes into
an array, in compressed format (making the matrix read-only from then
on).
But addressing elements using lists of indices seems a bit odd to me.
Is there some other way?
I have found some packages out there, but they seem to handle
2-dimensional matrices only.
Thank you,
J.
On Jun 1, 9:45 pm, Jeronimo Pellegrini <····@aleph0.info> wrote:
> Hello.
>
> I was wondering what would be an efficient way to implement
> three-dimensional sparse matrices in Common Lisp (from what I see in the
> HyperSpec, there is no native support for them).
> I also need them to have fast access...
>
[...]
> But addressing elements using lists of indices seems a bit odd to me.
> Is there some other way?
You can use recursive (eql) hash tables. So the first index would be
associated with a hash table, where you would then look the second
index up etc until the last one. The code is also a bit more complex
because you have to create the hash tables as needed.
Paul Khuong
From: Jeronimo Pellegrini
Subject: Re: Sparse matrices in Common Lisp? (with 3 dimentions)
Date:
Message-ID: <f3rclg$e6o$1@aioe.org>
On 2007-06-02, Paul Khuong <·······@gmail.com> wrote:
>> I was wondering what would be an efficient way to implement
>> three-dimensional sparse matrices in Common Lisp (from what I see in the
>> HyperSpec, there is no native support for them).
>> I also need them to have fast access...
>>
> [...]
>> But addressing elements using lists of indices seems a bit odd to me.
>> Is there some other way?
> You can use recursive (eql) hash tables. So the first index would be
> associated with a hash table, where you would then look the second
> index up etc until the last one. The code is also a bit more complex
> because you have to create the hash tables as needed.
Yes, but I'll have these inside deply nested loops, so I'd rather have
the best possible access time. Having to look up three hash tables is
much slower than oe or two direct accesses to an array, isn't it?
(I could be wrong...)
J.
On 2 Jun, 10:18, Jeronimo Pellegrini <····@aleph0.info> wrote:
> On 2007-06-02, Paul Khuong <·······@gmail.com> wrote:
>
> >> I was wondering what would be an efficient way to implement
> >> three-dimensional sparse matrices in Common Lisp (from what I see in the
> >> HyperSpec, there is no native support for them).
> >> I also need them to have fast access...
>
> > [...]
> >> But addressing elements using lists of indices seems a bit odd to me.
> >> Is there some other way?
> > You can use recursive (eql) hash tables. So the first index would be
> > associated with a hash table, where you would then look the second
> > index up etc until the last one. The code is also a bit more complex
> > because you have to create the hash tables as needed.
>
> Yes, but I'll have these inside deply nested loops, so I'd rather have
> the best possible access time. Having to look up three hash tables is
> much slower than oe or two direct accesses to an array, isn't it?
> (I could be wrong...)
>
> J.
Google ate my last response, so here's the short version.
If you can create a structure that places elements used together in
the same inner loop, together in the same inner most hash table, you
can expect significant speed ups due to caching. Working with
inference algorithms on Markov random fields, I managed to speed up a
single iteration from the best part of an hour to under 1 minute just
by doing this.
If you're working with uniquely structured sparse data sets and
concerned by speed, it's probably worth rolling you own accessors out
of the basic components.
(crosses fingers and hopes the mail gets through)
Chris
From: Mark Hoemmen
Subject: Re: Sparse matrices in Common Lisp? (with 3 dimentions)
Date:
Message-ID: <f3r8hi$23f2$1@geode.berkeley.edu>
Jeronimo Pellegrini wrote:
> I was wondering what would be an efficient way to implement
> three-dimensional sparse matrices in Common Lisp (from what I see in the
> HyperSpec, there is no native support for them).
> I also need them to have fast access...
When you say "three-dimensional sparse matrices," do you mean "sparse
3-D tensors"? I might be able to find someone for you who works on this
sort of thing, but I'll need some more information about the operations
you wish to perform with the tensors.
mfh
From: Jeronimo Pellegrini
Subject: Re: Sparse matrices in Common Lisp? (with 3 dimentions)
Date:
Message-ID: <f3rdrv$itk$1@aioe.org>
On 2007-06-02, Mark Hoemmen <············@gmail.com> wrote:
> Jeronimo Pellegrini wrote:
>> I was wondering what would be an efficient way to implement
>> three-dimensional sparse matrices in Common Lisp (from what I see in the
>> HyperSpec, there is no native support for them).
>> I also need them to have fast access...
>
> When you say "three-dimensional sparse matrices," do you mean "sparse
> 3-D tensors"? I might be able to find someone for you who works on this
> sort of thing, but I'll need some more information about the operations
> you wish to perform with the tensors.
That would be really nice -- thank you!
I need to store large matrices that represent some particular structure,
and run an algorithm that (non-destructively) produce something else.
The algorithm has matrix access inside nested loops, that's why I'd like
to have fast access.
The operations are very specific, and it's also where I'll do some
optimizations, so I'd be happy with the representation part only. :-)
(It's similar to the transition probability matrix in a very sparse
Markov chain, but it's 3D because there is one more parameter besides
previous state and next state.)
J.
From: Mark Hoemmen
Subject: Re: Sparse matrices in Common Lisp? (with 3 dimentions)
Date:
Message-ID: <f3t6i0$2l7q$1@geode.berkeley.edu>
Jeronimo Pellegrini wrote:
> On 2007-06-02, Mark Hoemmen <············@gmail.com> wrote:
>> When you say "three-dimensional sparse matrices," do you mean "sparse
>> 3-D tensors"? I might be able to find someone for you who works on this
>> sort of thing, but I'll need some more information about the operations
>> you wish to perform with the tensors.
>
> That would be really nice -- thank you!
>
> I need to store large matrices that represent some particular structure,
> and run an algorithm that (non-destructively) produce something else.
> The algorithm has matrix access inside nested loops, that's why I'd like
> to have fast access.
> The operations are very specific, and it's also where I'll do some
> optimizations, so I'd be happy with the representation part only. :-)
> (It's similar to the transition probability matrix in a very sparse
> Markov chain, but it's 3D because there is one more parameter besides
> previous state and next state.)
Our research group works on automatically tuned software for sparse
matrix-vector multiplication and similar operations: see e.g.,
http://bebop.cs.berkeley.edu/oski/
We target our data structures and algorithms for operations which
involve streaming through all the nonzero values of the sparse matrix,
several times, and performing some arithmetic operations involving each
nonzero entry of the matrix and the corresponding vector entry. Sparse
matrix-vector multiplication is the canonical example of this. None of
our algorithms modify the nonzero values of the matrix (i.e., we don't
do factorizations).
The runtime cost is divided somewhere between the bandwidth costs of
streaming the entries, and the latency costs of irregular accesses to
the vector(s). How it's divided depends on the nonzero pattern of the
matrix.
If you're doing something that looks like sparse matrix-vector
multiplication, I might be able to help you. But if there's no
arithmetic or vectors involved, our research probably isn't much help
for you.
Best,
mfh
From: Mark Hoemmen
Subject: Re: Sparse matrices in Common Lisp? (with 3 dimentions)
Date:
Message-ID: <f3t6l8$2l7r$1@geode.berkeley.edu>
Jeronimo Pellegrini wrote:
> On 2007-06-02, Mark Hoemmen <············@gmail.com> wrote:
>> When you say "three-dimensional sparse matrices," do you mean "sparse
>> 3-D tensors"? I might be able to find someone for you who works on this
>> sort of thing, but I'll need some more information about the operations
>> you wish to perform with the tensors.
>
> That would be really nice -- thank you!
>
> I need to store large matrices that represent some particular structure,
> and run an algorithm that (non-destructively) produce something else.
> The algorithm has matrix access inside nested loops, that's why I'd like
> to have fast access.
> The operations are very specific, and it's also where I'll do some
> optimizations, so I'd be happy with the representation part only. :-)
> (It's similar to the transition probability matrix in a very sparse
> Markov chain, but it's 3D because there is one more parameter besides
> previous state and next state.)
Our research group works on automatically tuned software for sparse
matrix-vector multiplication and similar operations: see e.g.,
http://bebop.cs.berkeley.edu/oski/
We target our data structures and algorithms for operations which
involve streaming through all the nonzero values of the sparse matrix,
several times, and performing some arithmetic operations involving each
nonzero entry of the matrix and the corresponding vector entry. Sparse
matrix-vector multiplication is the canonical example of this. None of
our algorithms modify the nonzero values of the matrix (i.e., we don't
do factorizations).
The runtime cost is divided somewhere between the bandwidth costs of
streaming the entries, and the latency costs of irregular accesses to
the vector(s). How it's divided depends on the nonzero pattern of the
matrix.
If you're doing something that looks like sparse matrix-vector
multiplication, I might be able to help you. But if there's no
arithmetic or vectors involved, our research probably isn't much help
for you.
Best,
mfh
From: Jeronimo Pellegrini
Subject: Re: Sparse matrices in Common Lisp? (with 3 dimentions)
Date:
Message-ID: <f3vt3g$2ce$1@aioe.org>
On 2007-06-03, Mark Hoemmen <············@gmail.com> wrote:
> Our research group works on automatically tuned software for sparse
> matrix-vector multiplication and similar operations: see e.g.,
>
> http://bebop.cs.berkeley.edu/oski/
>
> We target our data structures and algorithms for operations which
> involve streaming through all the nonzero values of the sparse matrix,
> several times, and performing some arithmetic operations involving each
> nonzero entry of the matrix and the corresponding vector entry. Sparse
> matrix-vector multiplication is the canonical example of this. None of
> our algorithms modify the nonzero values of the matrix (i.e., we don't
> do factorizations).
Well, that sounds like what I'd need. :-)
> The runtime cost is divided somewhere between the bandwidth costs of
> streaming the entries, and the latency costs of irregular accesses to
> the vector(s). How it's divided depends on the nonzero pattern of the
> matrix.
>
> If you're doing something that looks like sparse matrix-vector
> multiplication, I might be able to help you. But if there's no
> arithmetic or vectors involved, our research probably isn't much help
> for you.
Not exactly matrix-vector multiplication, but I do have to run
operations through all non-zeroes of vectors and matrices.
So... I went to that URL you sent, and there seems to be a C library
there, but no Lisp... :-) Is the Lisp stuff available also?
Thanks,
J.
From: Mark Hoemmen
Subject: Re: Sparse matrices in Common Lisp? (with 3 dimentions)
Date:
Message-ID: <f444jh$1ffn$1@geode.berkeley.edu>
Jeronimo Pellegrini wrote:
> Not exactly matrix-vector multiplication, but I do have to run
> operations through all non-zeroes of vectors and matrices.
I need a little more information in order to help you -- can you
describe the operation for me? (You're welcome to e-mail me if you
don't want to talk about it in public.) OSKI is for operations like
y := A * x,
y := A^T * x,
y := A^T * A * x,
in which A is a sparse matrix, and x and y are dense vectors.
> So... I went to that URL you sent, and there seems to be a C library
> there, but no Lisp... :-) Is the Lisp stuff available also?
OSKI is not written in Lisp because we have to make it available for new
architectures, for which a Lisp system may not be available. But you're
welcome to link it to Lisp via some combination of CFFI and SWIG. I
have done this myself for my Sparse Matrix Converter:
http://bebop.cs.berkeley.edu/smc/
mfh
From: Mark Hoemmen
Subject: Re: Sparse matrices in Common Lisp? (with 3 dimentions)
Date:
Message-ID: <f445uq$1fpi$1@geode.berkeley.edu>
Jeronimo Pellegrini wrote:
> Not exactly matrix-vector multiplication, but I do have to run
> operations through all non-zeroes of vectors and matrices.
I need a little more information in order to help you -- can you
describe the operation for me? (You're welcome to e-mail me if you
don't want to talk about it in public.) OSKI is for operations like
y := A * x,
y := A^T * x,
y := A^T * A * x,
in which A is a sparse matrix, and x and y are dense vectors.
> So... I went to that URL you sent, and there seems to be a C library
> there, but no Lisp... :-) Is the Lisp stuff available also?
OSKI is not written in Lisp because we have to make it available for new
architectures, for which a Lisp system may not be available. But you're
welcome to link it to Lisp via some combination of CFFI and SWIG. I
have done this myself for my Sparse Matrix Converter:
http://bebop.cs.berkeley.edu/smc/
mfh
From: Jeronimo Pellegrini
Subject: Re: Sparse matrices in Common Lisp? (with 3 dimentions)
Date:
Message-ID: <f464vb$pdd$1@aioe.org>
On 2007-06-05, Mark Hoemmen <············@gmail.com> wrote:
> Jeronimo Pellegrini wrote:
>> Not exactly matrix-vector multiplication, but I do have to run
>> operations through all non-zeroes of vectors and matrices.
>
> I need a little more information in order to help you -- can you
> describe the operation for me? (You're welcome to e-mail me if you
> don't want to talk about it in public.) OSKI is for operations like
>
> y := A * x,
> y := A^T * x,
> y := A^T * A * x,
Ah, I see.
I also need matrix-matrix multiplication, and access to individual
components of the matrices (sometimes in no specific order).
I think having a simple sparse representation (as discussed in the
thread) will probably be OK. I'll try to translate something from
C code.
>> So... I went to that URL you sent, and there seems to be a C library
>> there, but no Lisp... :-) Is the Lisp stuff available also?
> OSKI is not written in Lisp because we have to make it available for new
> architectures, for which a Lisp system may not be available. But you're
> welcome to link it to Lisp via some combination of CFFI and SWIG. I
> have done this myself for my Sparse Matrix Converter:
>
> http://bebop.cs.berkeley.edu/smc/
OK... But I'd like to avoid CFFI, as this would complicate portability
(would have to keep compiling/maintaining separate libraries besides the
Lisp code).
Thanks for the help!
J.
From: Mark Hoemmen
Subject: Re: Sparse matrices in Common Lisp? (with 3 dimentions)
Date:
Message-ID: <f484dp$2gv3$1@geode.berkeley.edu>
Jeronimo Pellegrini wrote:
> On 2007-06-05, Mark Hoemmen <············@gmail.com> wrote:
>> Jeronimo Pellegrini wrote:
>>> Not exactly matrix-vector multiplication, but I do have to run
>>> operations through all non-zeroes of vectors and matrices.
>> I need a little more information in order to help you -- can you
>> describe the operation for me? (You're welcome to e-mail me if you
>> don't want to talk about it in public.) OSKI is for operations like
>>
>> y := A * x,
>> y := A^T * x,
>> y := A^T * A * x,
>
> Ah, I see.
> I also need matrix-matrix multiplication, and access to individual
> components of the matrices (sometimes in no specific order).
> I think having a simple sparse representation (as discussed in the
> thread) will probably be OK. I'll try to translate something from
> C code.
I've got some sparse matrix-matrix multiplication code in C somewhere --
poke me if you're interested. Are you doing some kind of multigrid?
Then what you really want is a sparse triple product; I've got
specialized code for that, that doesn't create a temporary matrix.
You'll have to do the unfolding on it in order for it to work on more
general tensors.
mfh
From: Mark Hoemmen
Subject: Re: Sparse matrices in Common Lisp? (with 3 dimentions)
Date:
Message-ID: <f484p6$2gv5$1@geode.berkeley.edu>
Jeronimo Pellegrini wrote:
> On 2007-06-05, Mark Hoemmen <············@gmail.com> wrote:
>> Jeronimo Pellegrini wrote:
>>> Not exactly matrix-vector multiplication, but I do have to run
>>> operations through all non-zeroes of vectors and matrices.
>> I need a little more information in order to help you -- can you
>> describe the operation for me? (You're welcome to e-mail me if you
>> don't want to talk about it in public.) OSKI is for operations like
>>
>> y := A * x,
>> y := A^T * x,
>> y := A^T * A * x,
>
> Ah, I see.
> I also need matrix-matrix multiplication, and access to individual
> components of the matrices (sometimes in no specific order).
> I think having a simple sparse representation (as discussed in the
> thread) will probably be OK. I'll try to translate something from
> C code.
I've got some sparse matrix-matrix multiplication code in C somewhere --
poke me if you're interested. Are you doing some kind of multigrid?
Then what you really want is a sparse triple product; I've got
specialized code for that, that doesn't create a temporary matrix.
You'll have to do the unfolding on it in order for it to work on more
general tensors.
mfh
Jeronimo Pellegrini wrote:
> But addressing elements using lists of indices seems a bit odd to me.
> Is there some other way?
You could do what C does under the hood for dense arrays.
If the dimensions of the matrix are l, m, and n, and
the indices of an element are i, j, and k, you can
compute a numerical hash key (m*n)*i + n*j + k. It would
take up much less space than a list, and I think it would be
faster too, but it would be harder to change the dimensions
after you set them.
--
Dan
www.prairienet.org/~dsb/
From: Jeronimo Pellegrini
Subject: Re: Sparse matrices in Common Lisp? (with 3 dimentions)
Date:
Message-ID: <f3rcrj$e6o$2@aioe.org>
On 2007-06-02, Dan Bensen <··········@cyberspace.net> wrote:
> Jeronimo Pellegrini wrote:
>> But addressing elements using lists of indices seems a bit odd to me.
>> Is there some other way?
> You could do what C does under the hood for dense arrays.
> If the dimensions of the matrix are l, m, and n, and
> the indices of an element are i, j, and k, you can
> compute a numerical hash key (m*n)*i + n*j + k. It would
> take up much less space than a list, and I think it would be
> faster too, but it would be harder to change the dimensions
> after you set them.
Yes, that's what I thought at first -- but I'd then use a more flexible
approach (the one I described) and transform it into the C-like one
later, when the matrix doesn't need to be resized anymore.
J.
On Sat, 02 Jun 2007 01:48:39 -0500, Dan Bensen <··········@cyberspace.net> said:
| ...
| You could do what C does under the hood for dense arrays.
| If the dimensions of the matrix are l, m, and n, and
| the indices of an element are i, j, and k, you can
| compute a numerical hash key (m*n)*i + n*j + k. It would
But what is particularly special about C here (as opposed to Common
Lisp itself, or essentially even Fortran, for that matter)?
And wouldn't that be (i*m + j)*n + k?
---Vassil.
--
The truly good code is the obviously correct code.
> On Sat, 02 Jun 2007 01:48:39 -0500, Dan Bensen <··········@cyberspace.net> said:
> | You could do what C does under the hood for dense arrays.
> | If the dimensions of the matrix are l, m, and n, and
> | the indices of an element are i, j, and k, you can
> | compute a numerical hash key (m*n)*i + n*j + k. It would
Vassil Nikolov wrote:
> But what is particularly special about C here (as opposed to Common
> Lisp itself, or essentially even Fortran, for that matter)?
Maybe nothing, but I know it's done in C, whereas in Haskell, for
example, I don't think it's guaranteed, and many languages, e.g. OCaml
and most scripting languages, don't even support multidimensional
arrays. So it's not by any means universal. C is just an example
I happen to be familiar with.
> And wouldn't that be (i*m + j)*n + k?
Obviously that's a more efficient factoring, but it returns the
same value, and I thought completely separating the i, j, and k
terms might be a little clearer. So both expressions are correct,
and I think the separated form is useful for explaining the idea
to someone.
--
Dan
www.prairienet.org/~dsb/
On Sat, 02 Jun 2007 15:46:16 -0500, Dan Bensen wrote:
>> And wouldn't that be (i*m + j)*n + k?
>
> Obviously that's a more efficient factoring, but it returns the
> same value, and I thought completely separating the i, j, and k
> terms might be a little clearer. So both expressions are correct,
> and I think the separated form is useful for explaining the idea
> to someone.
Not only that, but a good compiler will probably un-factor it anyway,
because the factorization has just introduced an unnecessary serial
dependency on a temporary result. Particularly if any of n or m are
compile-time constants, which they must be if you're using the [][]
notation.
Cheers,
--
Andrew
On Sun, 03 Jun 2007 10:17:16 +1000, Andrew Reilly <···············@areilly.bpc-users.org> said:
| On Sat, 02 Jun 2007 15:46:16 -0500, Dan Bensen wrote:
|| [Vassil Nikolov wrote:]
||| And wouldn't that be (i*m + j)*n + k?
||
|| Obviously that's a more efficient factoring, but it returns the
|| same value, and I thought completely separating the i, j, and k
|| terms might be a little clearer. So both expressions are correct,
|| and I think the separated form is useful for explaining the idea
|| to someone.
| Not only that, but a good compiler will probably un-factor it anyway,
| because the factorization has just introduced an unnecessary serial
| dependency on a temporary result. Particularly if any of n or m are
| compile-time constants, which they must be if you're using the [][]
| notation.
Comparing (m*n)*i + n*j + k and (i*m + j)*n + k even when m and n
are compile-time constants, to be able to perform the two
multiplications in parallel is not necessarily that important in the
context of an actual program, since there will almost certainly be
other candidates for the ALU's attention at the same time [1]. On
the other hand, one frequently occurring optimization case is
multi-dimensional array access in nested loops, where the obvious
common subexpression elimination would make this factoring
consideration irrelevant anyway. (Assuming that the loops are
nested in the appropriate (row-major or column-major) way, of
course. I vaguely recall that some architectures have a
three-operand instruction that does one multiplication and one
addition, useful in such cases, but I do not recall if it would do
(X+Y)*Z or X*Y+Z.)
[1] an optimizing compiler for the suitable kind of architecture
will emit quite a hash of instructions when going for speed
Then, and more importantly, m and n need not in fact be compile-time
constants, which is frequently the case with library routines, for
example. Fortran has supported that for ages if not from the very
beginning, and Algol has that, too, if I am not mistaken (and, of
course, a few of the lisp languages). Even before C99, this could
be achieved in the C family with pointers to arrays, though somewhat
less efficiently, even if the [][] subscript notation in a somewhat
un-C-like way hides the increase in computational cost.
With regards to clarity, I believe that in fact the sequence i,
i*m + j, (i*m + j)*n + k, ((i*m + j)*n + k)*p + l, etc. is in fact
more illustrative of the way arrays of greater dimensions grow out
of arrays of lower dimensions.
Incidentally, good compilers need to be written by someone; a good
understanding of mathematics is a requirement for that job---the
right stuff, of course, not what grade school calls "mathematics",
though in all fairness Horner's scheme is taught in high school.
(What was the verdict on Mathematics vs. Good Programmer? Or is the
jury still out?)
---Vassil.
--
The truly good code is the obviously correct code.
On Jun 2, 9:17 pm, Andrew Reilly <···············@areilly.bpc-
users.org> wrote:
> On Sat, 02 Jun 2007 15:46:16 -0500, Dan Bensen wrote:
> >> And wouldn't that be (i*m + j)*n + k?
>
> > Obviously that's a more efficient factoring, but it returns the
> > same value, and I thought completely separating the i, j, and k
> > terms might be a little clearer. So both expressions are correct,
> > and I think the separated form is useful for explaining the idea
> > to someone.
>
> Not only that, but a good compiler will probably un-factor it anyway,
> because the factorization has just introduced an unnecessary serial
> dependency on a temporary result. Particularly if any of n or m are
> compile-time constants, which they must be if you're using the [][]
> notation.
>
> Cheers,
>
> --
> Andrew
The docucmentation of SparseLib for C++ describe in detail some
storage strategies for sparce matrix:
http://math.nist.gov/sparselib++/sparselib-userguide.pdf
From: Dimiter "malkia" Stanev
Subject: Re: Sparse matrices in Common Lisp? (with 3 dimentions)
Date:
Message-ID: <4661A25A.4080201@mac.com>
You could try using goedel numbers
(in-package "CL-USER")
(defun gn-hash (x y z)
(let ((dx (/ 1.0d0 (coerce x 'double-float)))
(dy (/ 1.0d0 (coerce y 'double-float)))
(dz (/ 1.0d0 (coerce z 'double-float))))
(* (expt 2 dx) (expt 3 dy) (expt 5 dz))))
(defun gn-test2 (ht x y z)
(setf (gethash (gn-hash x y z) ht)
(format nil "~A ~A ~A" x y z)))
(defun gn-test ()
(let ((ht (make-hash-table)))
(loop :for x :from 1 :to 10
do (loop :for y :from 1 :to 10
do (loop :for z :from 1 :to 10
do (gn-test2 ht x y z))))
(loop :for x :from 1 :to 10
do (loop :for y :from 1 :to 10
do (loop :for z :from 1 :to 10
do (format t "x=~A y=~A z=~A-> ~A~&" x y z
(gethash (gn-hash x y z) ht)))))))
#+nil
(gn-test)
Jeronimo Pellegrini wrote:
> Hello.
>
> I was wondering what would be an efficient way to implement
> three-dimensional sparse matrices in Common Lisp (from what I see in the
> HyperSpec, there is no native support for them).
> I also need them to have fast access...
>
> I thought I would implement them as hash-tables, with wrappers over
> these two:
>
> (setf (gethash (list i j k) my-matrix) value)
> (gethash (list i j k) my-matrix)
>
> And then have some "commit" function that will map the non-zeroes into
> an array, in compressed format (making the matrix read-only from then
> on).
>
> But addressing elements using lists of indices seems a bit odd to me.
> Is there some other way?
>
> I have found some packages out there, but they seem to handle
> 2-dimensional matrices only.
>
> Thank you,
> J.
>
(message (Hello 'Dimiter)
(you :wrote :on '(Sat, 02 Jun 2007 10:01:14 -0700))
(
DmS> You could try using goedel numbers
DmS> (defun gn-hash (x y z)
DmS> (let ((dx (/ 1.0d0 (coerce x 'double-float)))
DmS> (dy (/ 1.0d0 (coerce y 'double-float)))
DmS> (dz (/ 1.0d0 (coerce z 'double-float))))
DmS> (* (expt 2 dx) (expt 3 dy) (expt 5 dz))))
but how to find valid range for (x y z) where gn-hash is a perfect hash?
single-float version lacks precision quite soon:
CL-USER> (expt 2 (/ 1.0 65000))
1.0000106
CL-USER> (expt 2 (/ 1.0 65001))
1.0000106
this is very elegant, but i'm afraid it looses to integer hash both when
used as perfect hash, and as imperfect.
)
(With-best-regards '(Alex Mizrahi) :aka 'killer_storm)
"I am everything you want and I am everything you need")
On 2007-06-02 14:53:45 -0400, "Alex Mizrahi"
<········@users.sourceforge.net> said:
> but how to find valid range for (x y z) where gn-hash is a perfect hash?
(defun gn-hash-broken-p (x)
(= (gn-hash x x x) (gn-hash x x (1- x))))
CL-USER> (gn-hash-broken-p 114140001)
NIL
CL-USER> (gn-hash-broken-p 114140002)
T
On 2007-06-03 00:52:07 -0400, Raffael Cavallaro
<················@pas-d'espam-s'il-vous-plait-mac.com> said:
> On 2007-06-02 14:53:45 -0400, "Alex Mizrahi"
> <········@users.sourceforge.net> said:
>
>> but how to find valid range for (x y z) where gn-hash is a perfect hash?
>
> (defun gn-hash-broken-p (x)
> (= (gn-hash x x x) (gn-hash x x (1- x))))
>
> CL-USER> (gn-hash-broken-p 114140001)
> NIL
> CL-USER> (gn-hash-broken-p 114140002)
> T
oops
(defun find-gn-hash-limit (start step)
(loop for i from start by step
when (gn-hash-broken-p i)
return (if (= step 1) i
(find-gn-hash-limit (- i step) (/ step 10)))))
CL-USER 53 > (find-gn-hash-limit 85000000 1)
85128770
(message (Hello 'Raffael)
(you :wrote :on '(Sun, 3 Jun 2007 01:49:52 -0400))
(
RC> CL-USER 53 > (find-gn-hash-limit 85000000 1)
RC> 85128770
by the way, via combinatorial/information-theoretic approach, it's possible
to prove that for 64-bit double-floats coords with range 0..4194303 will be
broken.
indeed, 64-bit floats can encode not more than 2^64 = 18446744073709551616
distinct tuples.
but if each coord is encoded with 22-bit integer (upto 4194304), there are
2^66 of possible tuples, so 2^66 - 2^64 = 55340232221128654848 tuples will
be ambiguous.
moreover, certianly this encoding has significant overhead, so 22-bit is a
very optimistic upper limit.
i think most tuples will be encoded to a number close to 1, so they'll have
same exponent. thus, if IEEE 754 encoding is used, 11 bits of exponent in
number representation are not used, as well as sign bit, so we have
effectively only 52 bits.
this gives us only 18 bits for each coord -- upto 262144.
limits of such scale makes this issue looking more practical.
)
(With-best-regards '(Alex Mizrahi) :aka 'killer_storm)
"I am everything you want and I am everything you need")
(message (Hello 'Raffael)
(you :wrote :on '(Sun, 3 Jun 2007 00:52:07 -0400))
(
??>> but how to find valid range for (x y z) where gn-hash is a perfect
??>> hash?
RC> (defun gn-hash-broken-p (x)
RC> (= (gn-hash x x x) (gn-hash x x (1- x))))
RC> CL-USER> (gn-hash-broken-p 114140001)
RC> NIL
RC> CL-USER> (gn-hash-broken-p 114140002)
RC> T
i was thinking in that direction first, but i'm afraid it's not that simple.
i'm afraid there are indentical triples of more complex structure, e.g. (x
(1+ y) (1- z)) vs (x y z)
or maybe ((+ x 2) (1- y) (1- z)) vs (x y z).
and i'd expect that such clashes will come much sooner than gn-hash-broken-p
detects, since change in one coord is compensate with change in other coord,
so it's more likely that the product will be same.
unfortunately i don't know a simple way to analyze this, except doing a
full scan..
)
(With-best-regards '(Alex Mizrahi) :aka 'killer_storm)
"I am everything you want and I am everything you need")
On 2007-06-03 06:06:44 -0400, "Alex Mizrahi"
<········@users.sourceforge.net> said:
> unfortunately i don't know a simple way to analyze this, except doing a
> full scan..
right, and it may even vary from one implementation to another.
On 2007-06-03 06:06:44 -0400, "Alex Mizrahi"
<········@users.sourceforge.net> said:
> unfortunately i don't know a simple way to analyze this, except doing a
> full scan..
Actually, a simple list turns out to be faster and cons less than
gn-hash, and there is no risk of collisions. The timings below are in
sbcl. With a dimension-limit of most-positive-fixnum (or even far less
as you noted) there are guaranteed to be hash collisions for gn-hash.
(defun point-hash-test
(&key (times 10) (key-generator #'list) (table-test #'equal)
(print t) (dimension-limit most-positive-fixnum))
(time (let ((temp-table (make-hash-table :test table-test)))
(format t "Timing for ~a:~%" key-generator)
(loop repeat times do
(let ((x (random dimension-limit))
(y (random dimension-limit))
(z (random dimension-limit)))
(declare (fixnum x y z))
(setf (gethash (funcall key-generator x y z) temp-table)
(list x y z))
(if print
(format t
"The point ~a ~a ~a produces the key:
~a
at which key the list ~a is stored.~%~%"
x y z
(funcall key-generator x y z)
(gethash (funcall key-generator x y z)
temp-table))
(gethash (funcall key-generator x y z) temp-table)))))))
(defun speed-test (times)
(point-hash-test :times times :key-generator 'list :table-test 'equal
:print nil)
(point-hash-test :times times :key-generator 'gn-hash :table-test
'eql :print nil))
CL-USER> (speed-test 10000)
Timing for LIST:
Evaluation took:
0.026 seconds of real time
0.025878 seconds of user run time
1.95e-4 seconds of system run time
0 calls to %EVAL
0 page faults and
4,184,320 bytes consed.
Timing for GN-HASH:
Evaluation took:
0.073 seconds of real time
0.065698 seconds of user run time
0.00704 seconds of system run time
[Run times include 0.023 seconds GC run time.]
0 calls to %EVAL
0 page faults and
5,938,088 bytes consed.
NIL
CL-USER>
On 2007-06-05 15:59:06 -0400, Raffael Cavallaro
<················@pas-d'espam-s'il-vous-plait-mac.com> said:
> and there is no risk of collisions.
The situation is worse than I thought. Re-reading the spec, sxhash
returns a non-negative fixnum, so you'd need nested hash tables anyway.
(defun point-hash-test
(&key (times 10)
(print t) (dimension-limit most-positive-fixnum))
(flet ((present-p (key hash-table)
(multiple-value-bind (value anyone-home) (gethash key hash-table)
(declare (ignore value))
anyone-home)))
(time (let ((temp-table-x (make-hash-table)))
(loop repeat times do
(let ((x (random dimension-limit))
(y (random dimension-limit))
(z (random dimension-limit)))
(declare (fixnum x y z))
(unless (present-p x temp-table-x)
(setf (gethash x temp-table-x) (make-hash-table)))
(unless (present-p y (gethash x temp-table-x))
(setf (gethash y (gethash x temp-table-x))
(make-hash-table)))
(setf (gethash z (gethash y (gethash x
temp-table-x))) (list x y z))
(if print
(format t
"The point
~a ~a ~a
is stored in nested hash tables as the list ~a.~%~%"
x y z
(gethash z (gethash y (gethash x temp-table-x))))
(gethash z (gethash y (gethash x temp-table-x))))))))))
in sbcl this takes about the same time as gn-hash with a single hash
table, but can represent sparse arrays where each dimension is
most-positive-fixnum. It does, however, cons more than twice as much.
CL-USER> (point-hash-test :times 10000 :print nil)
Evaluation took:
0.089 seconds of real time
0.08812 seconds of user run time
0.001038 seconds of system run time
[Run times include 0.038 seconds GC run time.]
0 calls to %EVAL
0 page faults and
13,774,256 bytes consed.
NIL
There's some scheme code for manipulating sparse matrices in two
forms,
and implementing sparse matrix multiplication, at
<http://www.math.purdue.edu/~lucier/615/software/sparse-linear-
algebra.scm>
It's written in Scheme + Gambit extensions + the Meroon object system
(which
at least has generics/methods/instances/...). The surrounding
diretory has
more code that shows how the stuff in sparse-linear-algebra.scm is
used.
It's not what you're looking for, but it may give you some ideas. And
sparse
matrix-vector multiply is just as fast as in C for my applications,
because
it's memory-bandwith bound.
Brad
From: Mark Hoemmen
Subject: Re: Sparse matrices in Common Lisp? (with 3 dimentions)
Date:
Message-ID: <f4a00i$308a$1@geode.berkeley.edu>
······@math.purdue.edu wrote:
> It's not what you're looking for, but it may give you some ideas.
> And
> sparse
> matrix-vector multiply is just as fast as in C for my applications,
> because
> it's memory-bandwith bound.
Sort of... ;-) (It depends on the matrix structure -- highly irregular
accesses to the source vector may throw off the cache.)
I'd be interested in seeing some performance results for C vs.
Scheme/Lisp with sequential sparse matrix-vector multiplication. It
would be nice to use Lisp macros for generating Lisp code, rather than
<insert your favorite scripting language here> scripts to generate C
code ;-P
You can also save some bandwidth by improving the data structure. For
CSR-style sparse matrix data structures, this usually consists of
compressing or otherwise encoding the index stream. For example, if
your nonzeros are arranged in small dense blocks, all of the same size,
then you need only store one index for each block. You can also divide
your matrix into larger blocks, such that 16-bit integer indices are
sufficient for the indices in each block.
We do some of these in OSKI:
http://bebop.cs.berkeley.edu/oski/
and we'd like to do more. We'll welcome your suggestions :)
mfh
(message (Hello 'Raffael)
(you :wrote :on '(Wed, 6 Jun 2007 13:09:51 -0400))
(
??>> and there is no risk of collisions.
RC> The situation is worse than I thought. Re-reading the spec, sxhash
RC> returns a non-negative fixnum, so you'd need nested hash tables anyway.
why?
hash-tables are designed to work with imperfect hash functions, each bucket
can contain more than one key, and implementation will compare keys with
EQUAL when finds appropriate bucket. i think that's quite OK until count of
non-zero entries in matrix is less or comparable with most-positive-fixnum.
dimension limit, as well as dimension count can be arbitrary.
i think on 32-bit machines with nearly 32-bit most-positive-fixnum there
will be not enough addess space to allocate most-positive-fixnum doubles
anyway, so hasing should be fine.
if on 64-bit machines sxhash is also (almost) 64-bit, it's fine too.
)
(With-best-regards '(Alex Mizrahi) :aka 'killer_storm)
"I am everything you want and I am everything you need")
On 2007-06-09 04:31:10 -0400, "Alex Mizrahi"
<········@users.sourceforge.net> said:
> hash-tables are designed to work with imperfect hash functions, each bucket
> can contain more than one key, and implementation will compare keys with
> EQUAL when finds appropriate bucket. i think that's quite OK until count of
> non-zero entries in matrix is less or comparable with most-positive-fixnum.
> dimension limit, as well as dimension count can be arbitrary.
As you say, it depends on how many non-zero's we have, but for a 3d
array - which is what the OP asked about - if each dimension is
most-positive-fixnum, a density of (/ (expt most-positive-fixnum 2)) -
which is quite sparse indeed - already gives us most-positive fixnum
non-zero entries.
From: Mark Hoemmen
Subject: Re: Sparse matrices in Common Lisp? (with 3 dimentions)
Date:
Message-ID: <f3t5up$2kn3$1@geode.berkeley.edu>
Jeronimo Pellegrini wrote:
> I was wondering what would be an efficient way to implement
> three-dimensional sparse matrices in Common Lisp (from what I see in the
> HyperSpec, there is no native support for them).
> I also need them to have fast access...
When you say "three-dimensional sparse matrices," do you mean "sparse
3-D tensors"? I might be able to find someone for you who works on this
sort of thing, but I'll need some more information about the operations
you wish to perform with the tensors.
mfh
Jeronimo Pellegrini <···@aleph0.info> writes:
> Hello.
>
> I was wondering what would be an efficient way to implement
> three-dimensional sparse matrices in Common Lisp (from what I see in the
> HyperSpec, there is no native support for them).
> I also need them to have fast access...
My PDE solving library Femlisp does have some of this, although this part
is not very optimized for speed. See the file
#p"femlisp:src;algebra;sparse-tensor.lisp" in the source tree which you can
obtain at www.femlisp.org or via CVS from the Savannah server. Here, the
sparse tensors are defined recursively as sparse vectors of sparse vectors.
Nicolas
(message (Hello 'Jeronimo)
(you :wrote :on '(Sat, 2 Jun 2007 03:45:19 +0200 (CEST)))
(
JP> I was wondering what would be an efficient way to implement
JP> three-dimensional sparse matrices in Common Lisp (from what I see in
JP> the HyperSpec, there is no native support for them).
JP> I also need them to have fast access...
do you know that premature optimization is the root of all evil? :)
JP> I thought I would implement them as hash-tables, with wrappers over
JP> these two:
JP> (setf (gethash (list i j k) my-matrix) value)
JP> (gethash (list i j k) my-matrix)
looks like a quite sane solution. you can make it a bit shorter via dotted
list.
JP> And then have some "commit" function that will map the non-zeroes into
JP> an array, in compressed format (making the matrix read-only from then
JP> on).
most likely when you'll run this on the real data, you'll see which parts
are slow. if 99% of time is spent dealing with compressed format, there's
absolutely no sense to optimize HT part, right?
but if your HT will be too slow, you can just try several different options
and pick fastest for your data/CL implementation/hardware.
)
(With-best-regards '(Alex Mizrahi) :aka 'killer_storm)
"I am everything you want and I am everything you need")