Hi,
> For those who've not read my latest blog, I'm working on implementing
> a new structure for the matrices module of sympy.
> I faced many decision problems while doing so, like what should be
> named what and what should be where.
>
> Should matrix and linear algebra be separated ?
Probably yes in favour of a clean software architecture.
> I plan to put only essential functionalities in the main class. Other
> functionalities will be put in a separate file matrixtools.py or
> matrixutils.py.
That is probably good but the difficulties arise at the boundaries.
What is still a member function and what is not ...?
> so, once again the questions are,
> What do you expect from the matrix module ?
Hmm, matrices are a fairly general "structure" so there are many things
I would expect from a rather complete implementation. Maybe you should
also look at the Maxima manual, for example at chapter "49. diag" and
"67. linearalgebra"
> What should be the public API ?
I see two "parts" of a public API. Some operations are more like member
functions while others are more independent of a concrete matrix and
should be implemented as free functions. (This is also where we go
more into the direction of procedural code.)
> What functionality should be given special status to reside inside the
> class ?
Every matrix for example has a shape, so functions for retrieving the
shape definitely should go inside the class.
Then there are more things that I think of being intrinsic to
a given matrix instance. As long as we talk about square matrices
sucht things include trace, determinant, eigenvalues, eigenvectors,
characteristic and minimal polynomial and probably a few more.
Further we should be able to reason about the usual properties like
if the matrix is symmetric, hermitian, skew, diagonal ...
And maybe things like definiteness but I'm not sure if that should
better go into the assumtions system.
Also I would expect member functions for the usual transformations
like transpose (including hermitian transpose), invert etc, so I like
M.T(), M.invert(), ...
instead of
transpose(M), invert(M), ...
But this also depends a little on the question if the matrix M
is mutable here. For mutables, the member functions are better,
but for immutable objects maybe the independent functions are
more natural.
If we go to rectangular matrices, some of the above functions
are not defined, so the question is how we want to handle this?
And we also get "new" functions, f.e. singular values, pseudo inverses
and all that generalisations. These would then also be member functions
and reduce to the "usual" functions if the matrix is square.
(Or do we want to have also a method to compute the singular values
of a square matrix?)
To show a code snippet:
def singularvalues(self, ...):
if not self.is_square():
return singularvalues
else:
# Two choices:
return eigenvalues
# or
return NotImplementedError
and
def eigenvalues(self, ...):
if self.is_square():
return eigenvalues
else:
# Two choices:
return singularvalues
# or
return NotImplementedError
In contrast to the intrinsic properties, there are for example the
decompositions like LU, QR, Schur, Jordan, Polar, ... decomposition
which yield several matrices. So free functions are more appropriate
here:
L, U, P = ludec(M)
Q, R = qr(M)
but already here we hit a corner case, namely the eigendecomposition,
I said above that eigenvalues and vectors are more intrinsic properties.
But whats about:
P, Lambda = eigendec(M)
where M = P^-1 * Lambda * P
As we get two matrices here there can (and maybe really should) as well
be a free function for this. For example we get the eigenvalues as a diagonal
matrix here contrary to M.eigenvalues() where we would get a list/dict with
values and multiplicities.
I talked about intrinsic operations above, and another example case to consider
is
the computation of null spaces. This can be seen as important as f.e. the
determinant
and maybe resides on the same level. But then the question is, how do we
compute it?
Because solving linear systems is in my view not a member function, thus:
y = lu_solve(A,x)
is "better" than:
y = A.lu_solve(x)
And in this course, we would end at:
N = nullspace(M)
So there are many examples where it's not that obvious I think.
A similar example are matrix norms. I said I consider determinants
to be tightly coupled to a matrix instance. But whats about norms?
M.frobenius()
And what happens when we want to add more norms? Do we really
have to extend the matrix class code for this? Is this good design?
Maybe not ...
By the way, some of the above is also dependent on the ground fiels,
so if we deal with real, complex, ... matrices ...
> For example, jacobian does not belong to the class, it belongs to
> specialmatrices.py, and should be in the linalg. namespace as
> linalg.jacobian
Things like the jacobian and hessian are not directly related to the
matrices. (If we, at least for the moment, ignore more theoretical
concepts like operators etc) They just happen to be matrices,
but this is also not always true! (Jacobian of f:R^n->R is a vector,
or a n times 1 matrix, Hessian of f:R^m->R^n is a tensor object.)
So these functions should rather seamlessly integrate with diff, Derivative
and indexed objects / tensors. And they should play along grad, div, curl
laplacian etc.
To make a last example, what's the jacobian of a matrix of integers?
I would argue that this makes no sense unless you implicitely reinterpret
integers as constant functions ...
An simple implementation of the jacobian could look like:
def jacobian(f, (x,y,z,...)):
return Matrix( [[ diff(f_i, x_j) for i in ...] for j in ... ] )
Another important topic is iteration. I think that matrices really
should support iteration. And there are at least three cases which
are useful:
for row in matrix.iterrows():
...
for col in matrix.itercols():
...
for cell in matrix.iterentries():
...
Of course I'm aware that not all can be supported with
the same efficiency depending on the storage layout.
Maybe this is less important for symbolic matrices which
are not that big?
Or there is an efficient (maybe lazy) transpose operation,
and then if f.e. iterrows() is efficient, itercols() replaced
by M.T.iterrows() is somewhat efficient too? I don't know ...
For Spare matrices maybe we want to have:
for cell in M.iternonzeroentries():
...
In the same direction we have slicing. And I think that
slicing is *really* important, maybe even more important than
iteration. So any matrix should support slicing like:
M[a:b, c:d] yielding a continuous submatrix
M[(a,b), (c,d)] yielding a submatrix consisting of M[a,c], M[a,d],
M[b,c] and M [b,d]
For block matrices with arbitrary (and potentially different) block sizes
this should work on block level of course.
So, this are some of my thoughts about the topic. I don't know
if everything makes sense, especially with your code. (I have not
looked into it!) But anyway I hope this gives some useful input.
-- Raoul
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To post to this group, send email to sympy@googlegroups.com.
To unsubscribe from this group, send email to
sympy+unsubscr...@googlegroups.com.
For more options, visit this group at
http://groups.google.com/group/sympy?hl=en.
