GSoC 2013 Application Saurabh Jha: Linear Algebra Module

Proposal is at http://www.google-melange.com/gsoc/proposal/review/google/gsoc2013/saurabh_jha/11001

Personal Background

Name : Saurabh Jha

University : Guru Nanak Dev University, Gurdaspur, India

Current Standing : Third Year, Bachelor of Technology(B.Tech)

Major : Computer Science and Engineering

Contact :

email : [email protected]

IRC : Saurabh

github : SaurabhJha (https://github.com/SaurabhJha)

Programming Background

I have one and a half year experience with python and for the past few months, I am writing idiomatic python, though I won't claim I have absolutely mastered python language. I also have some amount of experience writing C, Java and Scheme/Lisp I think python is a wonderful language because it facilitates in rapid development and the bugs are somewhat easier to fix in comparison to compiled languages like C and Java. I think what python lacks is low level control but since we are building a high level application, we don't need that much control over system resources. The speed of any computer algebra system comes from the algorithm it implements, not the language. The most advanced standard library I have used in python is a url fetching library urllib2 while working on a project of a web crawler. I am now fairly comfortable with git and github. I can now focus on actual work rather than the idiosyncrasies of the git.

Contributions to Sympy

https://github.com/sympy/sympy/pull/1452(Merged) https://github.com/sympy/sympy/pull/2021(WIP) https://github.com/sympy/sympy/pull/2052(WIP) Project

This proposal is based on discussions with Aaron Meurer in IRC as well as some documents[1],[2],[3],[4]

Description of Problem:

Sympy is written in an interpreted language Python. Since it is interpreted, there is a huge overhead of using it, esepecially while calling the functions and using advanced runtime features. A typical algorithm may have hundreds of function calls. This is causing some speed issues with the matrix module today. One possible way to make Matrix module fast is to get rid of all the object oriented structures and using only stand alone functions as in imperative programming. The problem with this is approach is that the use of this module in Sympy or in interpreter will become a nightmare.

So, we have two choices, either get the speed with difficult to use interface or get better user interface but with speed issues as it is now. One way to get the best of both worlds is to divide up these problems into seperate layers, with one layer dealing with speed and perforamance and other layer the usability. We can have a bottom up architecture where the speed issues are handled by having fundamental operations and functionality operating in bottom levels, which are not meant to be seen by the user and then these operations are basis of more user friendly functionality.

Other problems associated with existing infrastructure are described in [5], [6]

The module will also have increase in speed because dedicated algorithms for different representations of matrix and different ground types will be implemented, although in top levels the user need not worry about the representation and domain of his matrices, the job of figuring that out is of the module.

My proposal for Google Summer of Code 2013 aims to rewrite linear algebra module to have the characteristics of both these tradeoffs. We can have performance as well as usability if we divide this problem into multiple layers, just like it was done in polys module[7].

My proposal therefore is to rewrite linear algebra module to be more like the polys module, where the matrices have ground types exactly similar to the polys module and have different representation based on the type of the matrix. The matrices will be more abstract they are now, that is, the manipulation of matrices are independent of the internals of the matrices, domains and representations as long as matrices on operations are themselves compatible with each other(like, in addition, number of rows and columns should be same and in multiplication, number of columns in first matrix should be equal to number of rows of second matrix).

Domains: Domains will be the basic data types in our module. Mathematically, domians are defined by the arithmetic operations and the properties of those arithmetic operations that we can have for the elements of that domain. Therefore, the domains are implemented as classes with methods of basic arithmetic and their properties for that particular type. Domains are in three categories--

Abstract Domains Concrete Domains Implementation Domains Abstract Domains

Abstract Domains are not supposed to be instantiated. They collect the common properties of other domains. Abstract Domains have the following classes

Domain(object)--Implements the basic most basic and common operations Ring(Domain)--Implements the basic properties of Rings. Field(Ring)--Implements the basic properties of Fields SimpleDomain(Domain), base class for simple domains like Integers, Rationals etc Concrete Domains

Concrete Domains implements the acutal implementations. They inherit from Abstract Domains. Concrete Domains have the follwing classes

FiniteFields(Field, SimpleDomain)-- Implement Finite/Galois Fields Integering(Ring, CharacteristicZero, SimpleDomain)-- Implements Integers PolynomialRing(Ring)-- Implements properties of Polynomials RationalField(Field, CharacteristicZero, SimpleDomain) Implements properties of Ratinals FractionField(Field, CompositeDomain) Implements Fraction Field ExpressionDomain(Field, CharacteristicZero, SimpleDomain) Implements Expresssions RealDomain(CharacteristicZero, SimpleDomain) Implements Real Numbers Implementation Domains

Implementation Domains implements interfaces with python's basic types, gympy and mpmath.It consists of the following classes

PythonFiniteField(FiniteField, simpledomain) Finite Field based on python's integers PythonIntegerRing(IntegerRing) Integer ring based on python's integers GympyRationalField(RationaField) Rational field based on GMPY mpq class mpmathRealDomain(RealDomain) Domain for real numbers based on mpmath mpf type. Architecture:

This module will be internally organised as a four-level architecure. The bottom two levels will provide the common functionality with focus on speed and performance. The upper two levels uses the functionality to bring up a common user interface to the user, which will hide all the complexities and should be easy enough to use. Each level will have it's own syntax and structure. On bottom two levels attention will be paid to speed and correctness of algorithms. On above two levels attention will be paid to useability. Special attention will be paid that no duplication of algorithms takes place.

Level 0: Representation of Matrices

Now we will look at the Matrix represeantation. Matrix representation is divided up to Level 0 and Level 1 of the above architecture. The raw representation of matrices will be in Level 0 and will be discussed first. There will be following Matrix Representations--

Dense Matrix Dense Matrix will be implemented as a tuple of a list of lists containing individual rows of matrices.

For example, if we want to represent the 3 by 3 matrix containing rows

(1, 7, 9); (10, 19, 25) and (5, 7, 45). We will use

Matrix(1, 7, 9], [10, 19, 2], [5, 7, 45)

[1 7 9]

[10 19 2]

[5 7 45]

Sprase Matrix:: Sprase Matrix will be implemented as a list of (index, element) tuples, where element itself will be a tuple (row, column).

For example, if we want to represent

[0, 0, 6, 0, 9, 0, 0]

[2, 0, 0, 7, 8, 0, 4]

[10, 0, 0, 0, 0, 0, 0]

[0, 0, 12, 0, 0, 0, 0]

we would write

Sparse([((0, 2), 6), ((0, 4), 9), ((1, 0), 0), ((2, 4), 7).....]])

[0, 0, 6, 0, 9, 0, 0]

[2, 0, 0, 7, 8, 0, 4]

[10, 0, 0, 0, 0, 0, 0]

[0, 0, 12, 0, 0, 0, 0]

We would have a function for that will return the representation a one dimensional representation so that rows and columns tuples are sorted. That will be useful while performing arithmetic operations on matrices, because the arithmetic operations will go through much fewer loops than if the row and column tuple are already sorted(explained in next section) Vector :: Since Vectors are basically one dimensional matrices, we will use a simple list

Vector([1, 3, 5, 9, 11])

[1]

[3]

[5]

[11]

Level 0: Functionality

This is the lowest level of functionality which is not meant to be seen by the user. All the operations will be done in imperative programming with explicit passing of arguments, as opposed to autamatic passing done in object oriented programming and no runtime features are to be used. Also, no object oriented programming features will be used.

These constraints, although making things fast, will make it harder to implement. I am expecting that this level is going to take a comparitively large amount of time.

The algorithms will be grouped according to the types of matrices and their domain or ground type. The types that I am going to define are--

Sprase Matrix Dense Matrix Vectors We can demonstrate the workings of Level 0 by having a look at the algorithm of addition in sprase matrix

Determine the number of rows and columns of the matrix and check if they are equal Sort the tuples containing the indices in increasing order From first matrix to last element check if the index tuple matches any of the other tuple if(yes): add the corresponding elements of two matrices and append the result to resultant matrix else: append the element of the first matrix to resultant matrix

Roughly the following functions will be implemented in this level (notice the words functions as opposed to methods)--

Arithmetic operations (additions, subtraction, multiplication) Rank of Matrix Echloen form Reduced row echloen form Adjoint of Matrix Inverse of Matrix Hessian of Matrix Augmented Matrix

Level 1

This is the second level and is also not supposed to be used by the user but only to be used internally. This level wraps the functions present in the Level 0 in Object Oriented Manner. More specifically, it will create three separate classes for each of the matrix representations and will contain methods that will use the functions implemented in level 0. This is also the level where domains start playing the role One of the parameters of these classes will be domain. Based on that domain different implementations of common operations will be done based on the domain that is passed as an argument. To illustrate, let's see a very crude snippet of code implementing add functions for different domain.

class SparseM(matrix, domain):

def add(self, other):

   if (self is of domain ZZ):

          add the above algorithm described in Level 0, with addition module of ZZ

   if (self is of domain Polys):

           use the above algorithm described in Level 0 functionality, with addition module of Polys.

Note that the object oriented features will produce some overhead because of some implicit argument passing but since the functions used were implemented in procedural style, there will be significantly less overhead.

Level 1 will also have a parser that can determine from a a given matrix it's suitable representation and it's ground type. Ideally, we would like to insulate the user as much as possible from underlying details and present the user with a common interface with whatever matrix he can come up with. A very crude outline may be as follows

if (matrix has m by n elements and i number of elements are zero for some i < m*n):

 it is a sparse matrix

elif(matrix has only one column)

 it is a vector

else

 it is a dense matrix

Level 2

This level is the first one that the user will interact with. This level will unite the matrix methods for different classes corresponding to different representations of matrices into a single matrix class. This level will use the parser to determine the type of matrix and it's domain and will use the appropriate class from level 1, which will in turn use appropriate functions of level 0. An operation from this level can be demonstrated again with the help of addition operation.

Matrix.add(M1, M2):

if (M1 and M2 are spare matrix):

use the sparse class from level 1

if (M1 and M2 have polys elements):

   use polys method of sparse

          use addition of polys function from level 0.

Level 3

This function is the uppermost interface to the user. This method will allow the user to use global functions as in

M1 + M2

as opposed to

Matrix.add(M1, M2). This will allow to have a more appealing user interface and will be similar to other contemporary Computer Algebra systems. Here the methods of level 2 are wrapped up to have a single function.

For example let we have the method inverse in level 2 of matrix class. It's usage will be as

Matrix.inverse(M1)

result

However, if we define a public function Inverse using the above method, the inverse of a matrix can be written as

Inverse(M1)

result

This level along with it's functions will be form the public API of the Matrix module.

Summary

The whole project can be explained in a nutshell as follows--

In L0, first have a set representations of matrices define and implement functions for different representations separately without any use of OOP concept. In L1, define classes for different representation of matrices and use the functions of L0 for to implement methods of function In L1, implement a parser that given the matrix, determines it's ground type and appropriate representation In L2, clump together all the classes of L1 into a single Matrix class with different methods for different representations and ground types. Representation and ground types will be determined by parser In L3, define external functions using the methods of L2 to get a more procedural oriented and clean interface to the user, similar to other Computer Algebra Systems

Timeline

Each bullet point is intended to be a pull request. Community Bonding Period : Reading Matrix literature, studying polys module, checking inconsistencies in plan. week 1: Implementing abstract domains week 2: Implementing concrete domains week 3: Implementing implementation domains week 4: Implement different matrix representations week 5: Implementing Level 0(for sparse matrix) week 6: Implementing Level 0(for dense matrix and vectors) week 7: Implementing Level 1(for integer, rational and fraction, galois, fraction, real and polynomial) week 8: Implementing Level 1(Implementing Parser for determining matrix representation and domain types) week 9: Implementing Level 2(Implementing unified representation of matrices and union of all the oprations of level 1) week 10: Implementing Level 3(Equivalent public functions for all the methods of level 2) week 11: Implementing Level 3 week 12: Fixing inconsistencies, writing documentation

References

http://groups.google.com/group/sympy/browse_thread/thread/bda446046424ea4a/f0923435d11793b8?lnk=gst&q=Linear+Algebra+polys# http://groups.google.com/group/sympy/browse_thread/thread/c62066049ac74716/c108d468055a2568? lnk=gst&q=Linear+Algebra+polys#c108d468055a2568 http://groups.google.com/group/sympy/browse_thread/thread/74b9ba9f22e5259f/6bafee5d1a969879?lnk=gst&q=Linear+Algebra+polys#6bafee5d1a969879 http://groups.google.com/group/sympy/browse_thread/thread/64197e1a1257a31e/e86ccb457a7e6388?lnk=gst&q=Linear+Algebra+polys#e86ccb457a7e6388 https://code.google.com/p/sympy/issues/detail?id=1850 https://github.com/sympy/sympy/pull/1840 http://mattpap.github.io/masters-thesis/html/src/internals.html