Fwd: Re: GSoC: Fast Linear Algebra over Extension Fields

[Martin Albrecht]

Hi,

lets keep this on lela-users so others can weight in as well. I'll get back to
you with comments (if I have any :)) later today.

Cheers,
Martin

---------- Forwarded Message ----------
Subject: Re: GSoC: Fast Linear Algebra over Extension Fields
Date: Tuesday 23 Apr 2013
From: Fonyó Dávid <fony...@gmail.com>
To: martinr...@googlemail.com


Hello Martin,

I think I've written the patch what you asked for GSoC, but I don't know
what to do now. Here is my code. I represented the two matrix over Z[p]/(f)
as you wrote, and I made a function that compute the product of the
matrices. You didn't write that I should multiply the matrices, but this
seemed the most reasonable, so I hope that I didn't misunderstand you. I
tested it on several inputs. It has no type checking, should I write it? Is
there something that I should correct or rewrite?
Thank you,

David

#include <vector>
#include <lela/integer.h>
#include <lela/ring/modular.h>
#include <lela/matrix/dense.h>

using namespace LELA;

template<class Ring>
void MyFunction(const Ring& R,
std::vector<DenseMatrix<typename Ring::Element> >& C, //m x l
const std::vector<DenseMatrix<typename Ring::Element> >& A, //m x n
const std::vector<DenseMatrix<typename Ring::Element> >& B, //n x l
const std::vector<typename Ring::Element>& P) //I supposed that P has
leading coefficient 1
{
std::vector<typename Ring::Element> temp(A.size()+B.size(), R.zero());
//it'll contain the entry of C before the factorization with P
for(unsigned int row=0; row < C[0].rowdim(); ++row) {
for(unsigned int col=0; col < C[0].coldim(); ++col)
{
for(unsigned int t=0; t < temp.size(); ++t)
{
R.copy(temp[t],R.zero());
}
for(unsigned int k=0; k < A[0].coldim(); ++k)
{
for(unsigned int i=0; i < A.size(); ++i) {
for(unsigned int j=0; j < B.size(); ++j) {
R.axpyin(temp[i+j],A[i][row][k],B[j][k][col]);
}
}
}
//factorize with p
unsigned int i = temp.size()-1;
while(i >= P.size()-1)
{
unsigned int j=i+1-P.size();
//std::cout<<"j"<<j<<std::endl;
typename Ring::Element c;
R.neg(c,temp[i]);
for(unsigned int k=0; k<P.size(); ++k)
{
R.axpyin(temp[j+k],P[k],c);
}
--i;
}
for(unsigned int i=0; i < P.size()-1; ++i)
{
R.copy(C[i][row][col],temp[i]);
}
}
}
}



2013/4/15 Martin Albrecht <martinr...@googlemail.com>

> Hey Dávid,
>
> (this should at least also be on lela-users)
>
> It would be very helpful if you could provide a little patch to LELA just
> to
> see whether you find your way round the library etc.
>
> For example, you could take two std::vector of mod p matrices and consider
> these vectors as matrices with polynomial entries. Then, you could perform
> schoolbook quadratic polynomial multiplication on these polynomials and
> take
> the result modulo some minimal polynomial represented as a std::vector<int>
> (or so).
>
> Something like that.
>
> Let me know if this is unclear and ask for help on lela-users if you get
> stuck.
>
> On Monday 15 Apr 2013, Dávid Fonyó wrote:
> > Hello,
> >
> > My name is Dávid Fonyó. I'm a third year Mathematics BSc and second year
> > Computer Science BSc student in Eötvös Loránd University, Budapest,
> > Hungary. I would like to join the Google Summer of Code 2013 program, and
> > I'm really interested in the "Fast Linear Algebra over Extension Fields"
> > project.
> >
> > As a prospective mathematician I've learned a lot of things that help me
> > understand the mathematical background of algorithms:
> > Algebra (4 semesters), Operation research (2 semester), Theory of
> > Computation, Numerical methods (2 semester), etc.
> > I have a good experience using C++, and I'm familiar with Java, Pascal,
> Ada
> > and Python.
> >
> > I've already checked the listed algorithms and now I'm reading the
> > references. I would like to ask what the next step will be? What kind of
> > patch should I have to write?
> > Thank you!
> >
> > Best regards,
> > Dávid Fonyó
>
> Cheers,
> Martin
>
> --
> name: Martin Albrecht
> _pgp: http://pgp.mit.edu:11371/pks/lookup?op=get&search=0x8EF0DC99
> _otr: 47F43D1A 5D68C36F 468BAEBA 640E8856 D7951CCF
> _www: http://martinralbrecht.wordpress.com/
> _jab: martinr...@jabber.ccc.de
>
> --
> --
> You received this message because you are subscribed to the Google
> Groups "lmnd-devel" group. Visit this group at
> http://groups.google.com/group/lmnd-devel?hl=en
> lmonade: http://www.lmona.de
>
> ---
> You received this message because you are subscribed to a topic in the
> Google Groups "lmnd-devel" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/lmnd-devel/Mm06g0qfNUQ/unsubscribe?hl=en
> .
> To unsubscribe from this group and all its topics, send an email to
> lmnd-devel+...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>
>



--
name: Martin Albrecht
_pgp: http://pgp.mit.edu:11371/pks/lookup?op=get&search=0x8EF0DC99
_otr: 47F43D1A 5D68C36F 468BAEBA 640E8856 D7951CCF
_www: http://martinralbrecht.wordpress.com/
_jab: martinr...@jabber.ccc.de

[Martin Albrecht]

Hi,

so, from what I can tell your solution is correct. However, it kind of misses the point, which means I did not explain this properly. Sorry about that!

You are doing the multiplications "by hand", i.e., you coded the usual triple-loop. However, the idea is to use the fast implementations of matrix multiplication available for GF(p). To illustrate this point, consider the following Sage code:

m,n,l,k = 2,2,2,2

A = [random_matrix(GF(7), m, n) for _ in range(k)]

B = [random_matrix(GF(7), n, l) for _ in range(k)]

C = [matrix(GF(7), m, l) for _ in range(2*k-1)]

for i in range(2):

    for j in range(2):

        C[i+j] += A[i]*B[j]

print "using fast arithmetic"

print sum([C[i]*x^i for i in range(2*k-1)])

print

P.<x> = GF(7)[]

print "using slow arithmetic"

Ax = sum([A[i]*x^i for i in range(k)])

Bx = sum([B[i]*x^i for i in range(k)])

Cx = Ax * Bx

print Cx

I put a copy here, so you can play with it:

http://aleph.sagemath.org/?z=eJx9j8tOwzAQRfeW_A9XXcWpQaIbJESRYkuwZW-FKoK2TBI7yHGFPx87UAR9jVcenzvHY6WTveywxEJOh7MqXYxv3NtgV7YJnmLx9FjcCgkr4QQ2g8cK5JCY7broRM2ZOptxEv3pjM6ZI8ExvCi7q5scyH360xd3nCFV7rcn-rm0oXlbY75EZagulWnToA9PLmC2G8ltsWnGgMZTeLfrQK-z_fO4s4XRORRfCIfu7z-JH5gzzp6v7-NDWmlaxRxaxn74_G-pYoInSXVG0k0C9cupi5zOXBpaQsW9XMcv57uMdw==&lang=sage 

Not sure whether it becomes clear but the first implementation uses the fast code available for matrices over GF(p), while the second uses the generic slow code for computing with matrices with polynomial coefficients.

I hope this helps, ping me on this list if you have further questions or I should explain this more.

[Dávid Fonyó]

Hi,

I rewrote it, I think this time it's what you asked. It makes the same result.

David

template<class Ring>

void Foo(const Ring& R,

std::vector<DenseMatrix<typename Ring::Element> >& C,

const std::vector<DenseMatrix<typename Ring::Element> >& A,

const std::vector<DenseMatrix<typename Ring::Element> >& B,

const std::vector<typename Ring::Element>& P)

{

std::vector<DenseMatrix<typename Ring::Element> > tempC (A.size()+B.size()-1,

DenseMatrix<typename Ring::Element> (A[0].rowdim(), B[0].coldim()));

Context<Ring> ctx(R);

for(unsigned int i=0;i<tempC.size(); ++i) {BLAS3::scal(ctx,R.zero(),tempC[i]);}

for(unsigned int i=0; i < A.size(); ++i) {

for(unsigned int j=0; j < B.size(); ++j) {

BLAS3::gemm(ctx,R.one(),A[i],B[j],R.one(),tempC[i+j]);

}

}

for(unsigned int row=0; row < C[0].rowdim(); ++row) {

for(unsigned int col=0; col < C[0].coldim(); ++col)

{

unsigned int i = tempC.size()-1;

while(i >= P.size()-1)

{

unsigned int j=i+1-P.size();

typename Ring::Element c;

R.neg(c,tempC[i][row][col]);

for(unsigned int k=0; k<P.size(); ++k)

{

R.axpyin(tempC[j+k][row][col],P[k],c);

}

--i;

}

}

}

C.resize(P.size()-1);

for(unsigned int i=0; i < C.size(); ++i)

{

BLAS3::copy(ctx,tempC[i],C[i]);

}

}

[Martin Albrecht]

Almost :) I'd suggest to use matrix additions to do the modular reduction in
the end instead of handcoding it as a double loop. But in any case, I see that
you got the point.

For the future: could you either link to a github or bitbucket repository or
attach .cpp files: copy'n'pasting code snippets + writing the boilerplate
around it, when you have perfectly functioning standalone programme is a bit
pointless.

I guess it's time to prepare that proposal? :)

> >> https://groups.google.com/d/topic/lmnd-devel/Mm06g0qfNUQ/unsubscribe?hl=

> >> en
> >>
> >> > .
> >> > To unsubscribe from this group and all its topics, send an email to
> >> > lmnd-devel+...@googlegroups.com.
> >> > For more options, visit https://groups.google.com/groups/opt_out.

[Dávid Fonyó]

Hi,

you're right, I rewrote the modular reduction part.

Actually, I'm new with github, but I created an account, and here is my code:

https://github.com/fonyodav/GSoC-2013/blob/master/lela_test

What should I do now?

[Martin Albrecht]

> What should I do now?

You'd now write a project proposal for the lmonade project (feel free to ask
for feedback here), deadline: May 3rd.

> Actually, I'm new with github, but I created an account, and here is my
> code:
> https://github.com/fonyodav/GSoC-2013/blob/master/lela_test

Note that this is not a in a usable form: this is neither a complete/stand-
alone C++ file nor does the file even have the right ending (.cpp). You must
have this code in a usable form because you tested it, why don't you just
commit that to your github repository?

On Thursday 25 Apr 2013, Dávid Fonyó wrote:
> Hi,

> you're right, I rewrote the modular reduction part.

[Dávid Fonyó]

Hi,

the LELA tutorial says that we can initialise matrices with a VectorStream. But I can't find a reference for that. Could you help me? Thank you!

David

[Bradford Hovinen]

Hello Dávid,

There are some examples in tests/test-blas-level3.h in the function testBLAS3ModulesConsistency. You can create a VectorStream over a particular field to generate up to a given number of random vectors of a given length and then feed that stream into the constructor of a matrix. The rows of the matrix are then automatically filled with the generated vectors.

Best regards,

Bradford Hovinen

2013/5/6 Dávid Fonyó <fony...@gmail.com>

--
You received this message because you are subscribed to the Google Groups "LELA Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lela-users+...@googlegroups.com.

[Dávid Fonyó]

Thank you Bradford!

Martin, I've revised my code, it's on github:

https://github.com/fonyodav/GSoC-2013

Dávid