Linear Algebra Vision

Introduction

Linear algebra is an important part of mathematics and a very important part of computational mathematics. The number of linear algebraic projects started within SymPy reflects this. Unfortunately this enthusiastic development has resulted in a number of very good but very separated projects. This article outlines a broad top-down vision of what we would like in a Linear Algebra system

Note: this article so far has only a single author. Help by adding your thoughts.

Needs

SymPy has the following needs for linear algebra

External

• Represent and manipulate small dense matrices with explicit entries (currently Matrix)
• Represent large sparse matrices with explicit entries (currently SparseMatrix?)
• Standard mathematical operations on these two types (i.e. solve, eigenvectors)
• Represent these matrices inside SymPy objects (ImmutableMatrix in PR 1015)
• Represent purely symbolic matrices, i.e. (X*Y).T (currently MatExpr)
• Matrix Calculus (i.e. d(xT A x)/ dx)
• Everything above, but now replace "Matrix" with "Tensor"
• Tensor index notation (What is a good summation notation easily expressible in Python?)
• Implement some well known tensor math - i.e. Relativity, fluids
• Interact with numpy - NumPy objects can contain SymPy ones. SymPy Matrices can not contain NumPy arrays
• Code generation (somewhat implemented in tensor.Indexed)
• Fancy indexing? How much of numpy functionality can we bring over?
• A common case in scientific computing is to have a linear operator rather than a matrix. I.e. an object for which individual indices are not available but the "multiply by a vector" operation is readily available.

Internal

Matrices are used internally by SymPy modules. These are currently using Dense Matrix objects but could benefit from a good efficient sparse solver

• Integration (indefinite)
• Solvers

Known issues

Mutability

For algorithmic performance some Matrix types will need to have rapidly changing entries. SymPy Basic objects must be immutable. This forces a split in the type system. Good design must be planned out so that this split can be executed with a minimum of code complexity.

Complex Interaction

Matrices and Tensors are used in a great variety of contexts. Building a simple system which can handle all of these contexts in a unified way may be challenging. Care must be taken to avoid code complexity.

Replication of Expr code

For pure symbolics (i.e. MatrixSymbol) we must build a system to represent and simplify tensor/matrix expressions. This system will be very similar but not identical to the Expr subclasses (Add, Mul, One, Zero, etc....) How can we make our own Add/Mul/Identity/Zero objects and robust simplification routines without redoing all of this work?

Proposals

To accomplish this while minimizing code complexity I think we need to look for clean places to conceptually break up this problem.

Storage vs. Algorithms

First, lets separate data storage from algorithms in explicitly represented matrices

Storage

We propose the following storage types

• DenseNDArray - SubContainer types
• --- list - Mutable
• --- tuple - Immutable
• --- numpy.ndarray - Fast, purely numeric. Use flags for immutability.
• SparseNDArray -
• --- dict - Use flags for immutability.
• --- scipy.sparse (fast but 2d matrices only)
• --- List-of-lists-of-lists-of...?
• Purely functional

Algorithms

Many algorithms can be performed without knowledge of the underlying storage. Merely having index access is sufficient. For a good list of these look into sympy/matrices/matrices.py

One of our goals though is for fast Sparse Matrices. For this it seems we may need to have a few special methods which are aware of underlying representation.

Symbolics

Tensor Indexing

Tensor calculus asks that we write equations like the following

$$R_{ab} - \frac{1}{2}g_{ab}R + g_{ab}\Lambda = \frac{8\pi G}{c^4}T_{ab}$$

We would like to be able to express this in the SymPy language easily and elegantly. What is the best way to get equations out of this out of the Python syntax?

What about contravariance/covariance? How should we represent different types of indices

One option

i,j,k,a,b,c,d = indicies('ijkabcd')
g[i,j]*u[0,i] # $g_i^j * u^i$
R[a|b, c|d] # $R_{ab}^{cd}$

Matrix should be a special case of Tensor. It would be nice if Matrix.matrix_multiply could be written as

def matrix_multiply(self, other):
    i,j,k = indicies('ijk')
    return self[i,j]*other[j,k]