GSoC 2022: Linear Algebra

[Oscar René Garzón Castro]

Hi all,

My name is Oscar Garzon. I have been programming in Python for two years, school projects and also personal projects. I amm a 3rd-year undergraduate student, studying mathematics at University of Mexico (UNAM). I've taken two courses in Linear algebra. I've used SymPy to solve equations systems.

I am interested in implementing Linear algebra and looking for a 350 hours project. For a project of this length, my proposal should include all the subtasks in Linear algebra or only in a number of them?

Kind Regards,
Oscar Garzon

[Aaron Meurer]

Can you clarify what you mean by this? SymPy already has a matrix
module, and there are several things that can be improved in it.

Aaron Meurer

>
> Kind Regards,
> Oscar Garzon
>
> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/2faf8c4e-f887-47ca-92a3-7d443ddad936n%40googlegroups.com.

[Alan Bromborsky]

If you are interested in linear algebra you might find the following book and link interesting -

https://smile.amazon.com/Linear-Geometric-Algebra-Alan-Macdonald/dp/1453854932/ref=sr_1_1?crid=386O85ZQLUMEZ&keywords=alan+macdonald+linear+algebra&qid=1649881349&sprefix=alan+mac+donald+linear+algebra%2Caps%2C35&sr=8-1

and

https://galgebra.readthedocs.io/en/latest/

--

[Oscar René Garzón Castro]

I mean the points in the Linear algebra subsection of Idea Prompts (https://github.com/sympy/sympy/wiki/GSoC-Ideas#idea-prompts). I am interested in some of this points but I'm not sure how many of them are expected for a 350 hours project.

[Oscar Benjamin]

On Thu, 14 Apr 2022 at 03:13, Oscar René Garzón Castro
<oscar....@ciencias.unam.mx> wrote:
>
> I mean the points in the Linear algebra subsection of Idea Prompts (https://github.com/sympy/sympy/wiki/GSoC-Ideas#idea-prompts). I am interested in some of this points but I'm not sure how many of them are expected for a 350 hours project.

The notes there are very much out of date. There is plenty of work to
be done along the "specialised data types" idea. That has already been
started as the DomainMatrix class which is now used internally by the
Matrix class. Most of the operations for Matrix don't make use of the
more efficient DomainMatrix methods though. Here's a simple example:

In [5]: M = randMatrix(10) # 10x10 matrix

In [6]: %time ok = M.inv()
CPU times: user 76 ms, sys: 0 ns, total: 76 ms
Wall time: 77.8 ms

In [7]: %time ok = M._rep.to_field().inv()
CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 2.46 ms

The difference is starker if you increase the size of the matrix e.g.
for 20x20 it's 30ms vs basically forever. For a matrix of only
integers like this the DomainMatrix method is much faster and should
always be used. Currently the Matrix class still uses its much slower
inv method though.

What needs to be done:

1. Make Matrix use more efficient DomainMatrix methods.
2. Add/complete/optimise DomainMatrix methods for important algorithms
(inv, rref, charpoly etc).
3. Explore and benchmark when different operations are better done
with some particular domain and consider alternative algorithms e.g.
when is it better to use an explicit inverse formula or row reduction.
4. Add helpers to efficiently convert to different domains and to
decide which domain should be used.
5. Make better use of block decompositions of matrices (i.e.
connected_components).

There's plenty to do there but you need to look at the code yourself
and try examples to understand what that is.

The fastest library I know of for these operations is flint which
should be made usable as an optional dependency through the
python_flint bindings.

In [1]: M = randMatrix(100)

In [2]: import flint

In [3]: fM = flint.fmpz_mat([[int(Mij) for Mij in Mi] for Mi in M.tolist()])

In [4]: %time ok = fM.inv()
CPU times: user 428 ms, sys: 8 ms, total: 436 ms
Wall time: 433 ms

In [5]: %time ok = M._rep.to_field().inv()
CPU times: user 6.97 s, sys: 4 ms, total: 6.98 s
Wall time: 6.98 s

You can see there that in this example flint is about 15x faster than
DomainMatrix. That's basically the difference between Python and C
when working with unbounded integers.

--
Oscar

> El miércoles, 13 de abril de 2022 a las 15:23:55 UTC-5, brombo escribió:
>>
>> If you are interested in linear algebra you might find the following book and link interesting -
>>
>> https://smile.amazon.com/Linear-Geometric-Algebra-Alan-Macdonald/dp/1453854932/ref=sr_1_1?crid=386O85ZQLUMEZ&keywords=alan+macdonald+linear+algebra&qid=1649881349&sprefix=alan+mac+donald+linear+algebra%2Caps%2C35&sr=8-1
>>
>> and
>>
>> https://galgebra.readthedocs.io/en/latest/
>>
>> On 4/12/22 8:43 PM, Oscar René Garzón Castro wrote:
>>
>> Hi all,
>>
>> My name is Oscar Garzon. I have been programming in Python for two years, school projects and also personal projects. I amm a 3rd-year undergraduate student, studying mathematics at University of Mexico (UNAM). I've taken two courses in Linear algebra. I've used SymPy to solve equations systems.
>>
>> I am interested in implementing Linear algebra and looking for a 350 hours project. For a project of this length, my proposal should include all the subtasks in Linear algebra or only in a number of them?
>>
>> Kind Regards,
>> Oscar Garzon
>>
>> --
>>
>> You received this message because you are subscribed to the Google Groups "sympy" group.
>> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
>> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/2faf8c4e-f887-47ca-92a3-7d443ddad936n%40googlegroups.com.
>
> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/5a62fe76-a526-43d8-bac1-cef6cd98835cn%40googlegroups.com.

[David Bailey]

Dear Group,

I would like to understand the workings of SymPy a little better, and to
that end, I thought it would be instructive to create a trivial 'special
function' that had some or all of the properties of the SymPy in-built
functions. The following link was helpful:

https://github.com/Shekharrajak/sympy/wiki/About-implementing-special-functions

However, I could not tell from that whether I could have the new
function stored separately from the directory tree holding the SymPy
software. I'd much prefer it if I could experiment without disturbing
the SymPy software at all.

I decided to call the function bailey, so following the recipe in the
above link, I created a file bailey.py in a directory outside the SymPy
tree, and to test it my calling Python from that directory, as follows:

from sympy.core.function import Function
class bailey(Function):
   """
   The experimental function
   """
   nargs=1
   @classmethod
   def eval(cls, arg):
           # This body does nothing useful at all
        z=arg[0]
        if z is S.Zero:
            return S.Zero
        else: return 1

Now I tried to use this module:

>>> import sympy
>>> import bailey

>>> bailey(a+3)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: 'module' object is not callable

Alternatively I tried:

>>> (a+3).bailey(1)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: 'Add' object has no attribute 'bailey'

I admit to being fairly confused at this level! What I would like to
know is:

1)       Can I experiment in this way without actually adding files to
SymPy itself or modifying anything there already?

2)       I would expect to get back a function call: bailey(a+3) or
perhaps the answer 1.

3)       Is there a place where this is described more completely?

Thanks,

David

[Aaron Meurer]

On Fri, Apr 15, 2022 at 2:23 PM David Bailey <da...@dbailey.co.uk> wrote:
>
> Dear Group,
>
> I would like to understand the workings of SymPy a little better, and to
> that end, I thought it would be instructive to create a trivial 'special
> function' that had some or all of the properties of the SymPy in-built
> functions. The following link was helpful:
>
> https://github.com/Shekharrajak/sympy/wiki/About-implementing-special-functions

This is a link to someone's fork of the wiki. The actual page is here
https://github.com/sympy/sympy/wiki/Development-Tips-:-About-implementing-special-functions

>
> However, I could not tell from that whether I could have the new
> function stored separately from the directory tree holding the SymPy
> software. I'd much prefer it if I could experiment without disturbing
> the SymPy software at all.

Yes, you can have it separate. The location of a function in the code
has no bearing on its behavior.

>
> I decided to call the function bailey, so following the recipe in the
> above link, I created a file bailey.py in a directory outside the SymPy
> tree, and to test it my calling Python from that directory, as follows:
>
> from sympy.core.function import Function
> class bailey(Function):
> """
> The experimental function
> """
> nargs=1
> @classmethod
> def eval(cls, arg):
> # This body does nothing useful at all
> z=arg[0]
> if z is S.Zero:
> return S.Zero
> else: return 1
>
> Now I tried to use this module:
>
> >>> import sympy
> >>> import bailey
>
> >>> bailey(a+3)
> Traceback (most recent call last):
> File "<stdin>", line 1, in <module>
> TypeError: 'module' object is not callable

When you run 'import bailey' then 'bailey' will be the module (because
you called your file bailey.py). To get the function do

from bailey import bailey

>
> Alternatively I tried:
>
> >>> (a+3).bailey(1)
> Traceback (most recent call last):
> File "<stdin>", line 1, in <module>
> AttributeError: 'Add' object has no attribute 'bailey'

You can't add new methods to existing SymPy objects. The first way you
tried is the correct way to access your bailey function.

>
> I admit to being fairly confused at this level! What I would like to
> know is:
>
> 1) Can I experiment in this way without actually adding files to
> SymPy itself or modifying anything there already?
>
> 2) I would expect to get back a function call: bailey(a+3) or
> perhaps the answer 1.

The function as you've written it will just return 1. If you want a
symbolic answer, eval() needs to return None (which is also the same
as not returning anything).

Since your eval() always returns something, the function will never be
symbolic. It is effectively no different from a normal Python function

def bailey(x):
if x == 0:
return 0
else:
return 1

>
> 3) Is there a place where this is described more completely?

I think the best docs right now are at
https://docs.sympy.org/latest/modules/core.html#id55. I will probably
end up writing something better at some point as this an important
thing that is somewhat underdocumented right now.

Aaron Meurer

>
> Thanks,
>
> David

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

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/04eaf831-b1a1-6f1c-68b7-2ae9a7d60c2e%40dbailey.co.uk.

[David Bailey]

Thanks for such an amazingly fast and detailed reply!

On 15/04/2022 21:32, Aaron Meurer wrote:
>
> Yes, you can have it separate. The location of a function in the code
> has no bearing on its behavior.

That is good!

> When you run 'import bailey' then 'bailey' will be the module (because
> you called your file bailey.py). To get the function do
>
> from bailey import bailey

I tried that too, but I got:

>>> bailey(a+3)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>

  File "C:\PythonSystem\lib\site-packages\sympy\core\cache.py", line
70, in wrapper
    retval = cfunc(*args, **kwargs)
  File "C:\PythonSystem\lib\site-packages\sympy\core\function.py", line
476, in __new__
    result = super().__new__(cls, *args, **options)
  File "C:\PythonSystem\lib\site-packages\sympy\core\cache.py", line
70, in wrapper
    retval = cfunc(*args, **kwargs)
  File "C:\PythonSystem\lib\site-packages\sympy\core\function.py", line
288, in __new__
    evaluated = cls.eval(*args)
  File "C:\SymPyProject\bailey.py", line 9, in eval
    # Value at zero
TypeError: 'Add' object is not subscriptable

>
>> Alternatively I tried:
>>
>> >>> (a+3).bailey(1)
>> Traceback (most recent call last):
>> File "<stdin>", line 1, in <module>
>> AttributeError: 'Add' object has no attribute 'bailey'
> You can't add new methods to existing SymPy objects. The first way you
> tried is the correct way to access your bailey function.
>
>> I admit to being fairly confused at this level! What I would like to
>> know is:
>>
>> 1) Can I experiment in this way without actually adding files to
>> SymPy itself or modifying anything there already?
>>
>> 2) I would expect to get back a function call: bailey(a+3) or
>> perhaps the answer 1.
> The function as you've written it will just return 1. If you want a
> symbolic answer, eval() needs to return None (which is also the same
> as not returning anything).
>
> Since your eval() always returns something, the function will never be
> symbolic. It is effectively no different from a normal Python function
>
> def bailey(x):
> if x == 0:
> return 0
> else:
> return 1

To be absolutely clear, the actual content of the function is just junk
at this point - I am just trying to get the basic mechanism to work.

David

[Aaron Meurer]

You are calling arg[0], but arg is already the argument of the
function. The [0] is typically only used when you have *args to accept
an arbitrary number of arguments and you want to get the first one.

Aaron Meurer

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

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/2e6860c8-b1dd-fc37-9fa5-5bfa4d74b8d3%40dbailey.co.uk.

[David Bailey]

Thanks again for a near instantaneous reply!

On 15/04/2022 23:39, Aaron Meurer wrote:
>
> You are calling arg[0], but arg is already the argument of the
> function. The [0] is typically only used when you have *args to accept
> an arbitrary number of arguments and you want to get the first one.

OK - that has fixed it after I corrected another minor error!

Just one possible use for this mechanism would be to  create a function
gamma() (say) which would contain all the original code from gamma in
SymPy, and then add calls to print to clarify how it works. I assume
that the line

import gamma from gamma

would replace the official version with this modification and everything
would work as expected?

David

[Aaron Meurer]

Yes, it ought to.

Aaron Meurer

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

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/b3d0f8e3-95ca-26e1-1021-9e91995fd883%40dbailey.co.uk.