Hi, I'm a little late to this conversation because I didn't follow it during the SciPy conference.
Comments inline: On Fri, Jul 15, 2011 at 10:01 AM, someone <[email protected]> wrote: > 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 ? So, here you need to be clear on what you think of as "linear algebra". I know what the difference is mathematically, but I am not entirely clear that you are meaning the same thing. To me, linear algebra is a theoretical way of looking at things. In the finite dimensional case, linear operators are equivalent to matrices. There are two main reasons why you would want to do "linear algebra" without directly using matrices: - Infinite dimensional vector spaces. This is probably not really much of a consideration in the CAS view, but to be sure, matrix representations only work in finite dimensional spaces. - Looking at things with respect to matrices requires you to have a basis. This is where I think some stuff that does direct things with linear operators without the user having to manually use matrices would be useful. For example, suppose we had a linear operator F(p) = p' + 2*p on the set Q[x] of polynomials in x. This set is infinite dimensional, but we can restrict ourselves to only polynomials of a certain degree (depending on the problem), so matrix representations are useful. But actually constructing the matrix of this operator is not trivial. Fist, you need a basis. Polynomials have a standard basis, but in some cases you might want to work in some other basis. Second, you must rewrite an arbitrary input in terms of this basis. Third, you must convert this into a Matrix. It would be awesome if you could just do something tell SymPy that that is your operator and it operates over the set, say P2 of polynomials of degree 2 or lower. This sort of thing might be beyond the scope of your project, but I just want to throw it out there. > > 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. Where things are is superficial, and can easily be changed without affecting the public interface. > > 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.) I also see two parts of the public API, but in a different sense. There is the user level stuff. This is pretty straight forward, as we already basically have something that works, and we just need to adapt the new stuff to work transparently with it. There is also the API for someone who wants to extend the module, say by writing a new matrix algorithm or representation. This is much less easy to get right. > >> 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. Mutability is an interesting point to raise here, because currently the mutable Matrix class is not changed by M.T or M.invert(). So actually the question arises as to exactly which things should be done mutably. Personally, I think something should be mutable only if it saves time or resources. If you have to basically reconstruct the Matrix, the operation might as well be immutable, as that is what people will expect anyway. I think this will always hold for inverse and for transpose it depends on the internal data type. Actually, that is something else to consider, that you cannot completely ignore the internal data type, as some operations are inherently more efficient in one way for one data type than in a completely different way for a different data type. The classic example of this is when you have an efficiency gain by doing things mutably. This is, e.g., why list.sort() modifies the list in place. > > 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 I would say raise NotImplementedError. Explicit is definitely better than implicit in this case. If singularvalues() return eigenvalues and eigenvalues() returns singular values, then it will be impossible for the user to tell which she is getting. > > 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. There is another side to this, which is that functions would pollute the namespace. Having functions in the namespace is not an inherently bad thing. The main problem is that these functions operate on Matrix objects, whereas most functions in SymPy operate on Expr objects. It depends on the name of course, but somethings like *_solve() and variants would be very confusing to the user against the normal SymPy solve() function. Of course, we could name things matrix_solve(), but at that point, you might as well just make them methods of Matrix. So I actually vote to have almost all if not all matrix operations as methods to Matrix. This might get hairy if there is some operation that has more than two operands, in which case just having a function might make sense. > > > 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(): > ... Should Matrix.__iter__ just iterate over the entries, so "for cell in matrix" is equivalent to "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 ... I think this depends on the internal representation of the Matrix. > > 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 > Aaron Meurer -- You received this message because you are subscribed to the Google Groups "sympy" group. To post to this group, send email to [email protected]. To unsubscribe from this group, send email to [email protected]. For more options, visit this group at http://groups.google.com/group/sympy?hl=en.
