Second draft of linearalgebra.mac for the contrib folder. #1185
Conversation
Update 06/12/23
Version 0.2 of linearalgebra.mac
Remove the s_test_case for vec_convert that seems to be causing issues.
The convenience functions `c` and `r` throw errors regularly if non-conforming vectors are added/multiplied inappropriately, so the original solution of directly converting to matrix form was inappropriate. The solution is to keep `c` and `r` as inert functions and use `vec_convert` to remove these as needed. Correspondingly, `vec_convert` now uses `errcatch` and will return the original expression when errors arise. A new predicate function, `vec_convertedp`, will detect whether an expression was successfully converted to matrix form.
Whilst editing questions I have found that saving and loading has been very burdensome compared to normal and my theory is that it runs all tests on every save, and on every preview (perhaps with deployed variants this would stop). This file exists with no `s_test_case` lines to see if the load times are improved.
Declared `c` and `r` to be nonscalar so that they behave appropriately with other matrices in expressions.
Fixed the vector display functions; they now display the tex-formatted entries rather than just the plaintex (how did I miss this earlier??)
Prevented divide by zero errors in integerify
Fixed subspace_equiv Added col_vecp and row_vecp
v0.2 of linearalgebra.mac features most of the functions I expect to use, and test cases for almost all of these functions. It contains: - Convenience functions for student vector input - Various predicate functions for vector and matrix properties - Functions to manipulate/extract parts of matrices - Functions to convert collections of vectors into a locally standard form (currently preferring lists of lists or matrices) - Functions to compare objects, typically checking for independence, or equivalence to subspaces. - Functions to convert collection of vectors to a basis - Extension of significantfigures to work for matrices - Extension of norm and cond functions to use the 2-norm - Extension of linsolve to accept matrix input and optionally use least squares - A function to produce smallest integer vector parallel to given vector (useful for eigenproblems!) - rowspace and nullTspace functions to pair with columnspace and nullspace - Rayleigh quotient calculation - Algebraic and geometric multiplicities of eigenvalues - Computing various matrices and factorisations, including rectangular diagonal matrix, projection matrices, QR, Jordan form, (orthgonal) diagonalisation, SVD, and pseudoinverse. Still to-do: - Possibly write some eigenvalue/eigenvector predicates (need to determine what is actually useful) - Work on the disp_eqns function that produces nice LaTeX for displaying systems of linear equations. - More functions will likely get added (and possibly removed) as work continues on local linear algebra course.
Copy of linearalgebra.mac v0.2.3 (committed 8th May 2024) with the s_test_case lines removed.
Implemented old fix for integerify (missed in the upgrade to 0.2.3).
Implemented old fix to integerify (missed in upgrade to 0.2.3)
There was a problem hiding this comment.
Ok, the fact that you (@LukeLongworth) have a _no_test.mac file vindicates my proposal in issue #1184 to move all the tests into a _test.mac file!
|
@LukeLongworth, thanks so much for pushing on with all of this. It looks super useful! I do like the We are about to refactor a lot of the Maxima code to include docs within the code, and then automatically generate the user docs from the docs in the code. That puts all the docs into one place. We might, at that point, move the test cases back into the main codebase but strip them out when we include the code within the question. I do really like the principle of having docs, code and tests all adjacent in one file. Perhaps we could have an email discussion about where you are with this? Would you like to include this in the next release and is (part?) of it ready? |
Following the fix to maths#1148: maths@a35ebc0
Deleting file with plan to rename `linearalgebra_no_test` to `linearalgebra`
Mostly a documentation update. Also added checks in lin_indp and subspace_equiv so that they don't throw errors when given lists of lists of unequal length.
| * @param[matrix] M a matrix | ||
| * @return[boolean] Does M have orthogonal columns? | ||
| */ | ||
| orthogonal_columnsp(M):= ev(diagp(transpose(M).M),simp); |
There was a problem hiding this comment.
Shouldn't the diagp be outside the ev?
I also noticed that sometimes, even with simp, complicated expressions (with square roots as they may appear in the context of orthonormal matrices) sometimes don't get fully simplified. For example,
only collapses into the identity matrix when applying expand. I wonder if it's worthwhile adding this?
There was a problem hiding this comment.
I think that the ev is fine in either position. I got a bit paranoid about simplification when writing, as I kept finding that certain functions or expressions were returning unevaluated unexpectedly. My interpretation of this particular expression is "take the expression diagp(transpose(M).M) and evaluate it with simplification turned on", so the transpose(M).M should still be evaluated as expected before the diagp predicate runs. My very quick check now suggests that both are fine, but if I see a counter-example I'm very happy to swap them.
Yes, that's a good point on expand. Could you please send me the problem case you've screenshotted above? I'll need to figure out which functions will need that and that looks like a good test case to use.
Thanks for your help with little project Georg! It's great to have a second pair of eyes on this :)
There was a problem hiding this comment.
Fair point about the ev, my intuition was that diagp returns a boolean anyway so I thought you'd want to simplify the expression within diagp, but I guess technically it might as well not evaluate the if-clauses and return them as is, for some reason, and if simplification is on, it'll apply recursively everywhere so no harm putting it on the outside? You probably have a better idea of how all of this works than I do...
The use case where I encountered this simplification issue is actually slightly more complex as I'm using a custom scalar product. I've been trying to reproduce the issue for the standard scalar product, but so far without success.
Anyway, below is some maxima code that you can run to see that expand sometimes does make a difference.
/* Compute orthogonal projection with respect to a custom scalar product */
proj(u, v, sp) := sp(v,u)/sp(u,u)*u;
/* Gram-Schmidt orthogonalization with respect to a custom scalar product */
gram_schmidt(M, sp) := block([u_1, u_2,u_3,v_2,v_3,AA],
/*Prepare Gram-Schmidt process*/
u_1: col(M,1),
v_2: col(M,2),
v_3: col(M,3),
/*orthogonalize*/
u_2: v_2 - proj(u_1, v_2, sp),
u_3: v_3 - proj(u_1, v_3, sp) - proj(u_2, v_3, sp),
/*normalize*/
u_1: u_1/sqrt(sp(u_1,u_1)),
u_2: u_2/sqrt(sp(u_2,u_2)),
u_3: u_3/sqrt(sp(u_3,u_3)),
AA: addcol(u_1, u_2, u_3),
return(AA)
);
/* A symmetric positive definite matrix and corresponding bilinear form (scalar product) */
A: matrix([3,-1,0],[-1,3,0],[0,0,2]);
bform(w,v) := transpose(w).A.v;
/* Orthogonalize the standard basis with respect to this scalar product */
baseM: ident(3);
OB: expand(gram_schmidt(baseM, bform));
/* Customized predicate functions */
orthogonal_columnsp(M, spMat):= diagp(ev(transpose(M).spMat.M,simp));
orthonormal_columnsp(M, spMat):= if matrixp(M) then is(ev(transpose(M).spMat.M,simp) = ident(second(matrix_size(M)))) else false;
/* Check outputs */
orthogonal_columnsp(OB, A);
orthonormal_columnsp(OB, A);
ev(transpose(OB).A.OB, simp);
ev(expand(transpose(OB).ev(expand(A.OB), simp)), simp);
- Added support for displaying augmented matrices easily using an inert function and texput. Also added `aug` and `de_aug` to convert matrices in and out of this mode. - Reworked `triu`, `tril`, `get_diag`, `triup`, `trilp`, `diagp`, and `diag_entries`. These should be significantly more efficient now. Thanks for the `genmatrix` tip-off Georg Osang! - Added `un_vec_convert` to change matrices into `c` and `r` form. Mostly intended for making teachers' answers the format we expect students to use. - Added a line to `make_list_of_lists` that will take a single vector and make it a list of a list. e.g. `[1,2,k]` -> `[[1,2,k]]`. - Edited `orthogonal_columnsp` and associated functions to forcibly expand the matrix before comparing. Thanks Georg! - Added support for eigenproblems. This includes - `eigenvaluep` and `eigenvectorp` which check whether a given value or vector is "eigen". Optionally checks them against a given vector or value respectively. - `get_eigenvalue` and `get_eigenvector` to extract the corresponding one from a given. Much nicer than using `eigenvectors` when you only want one! Can also orthonormalise the basis if requested. - Moved `Rayleigh` up to be with these. - Added support for vector parametric equations. It's a bit sparse, but these functions can check whether a given expression is of the desired form, and if so extract the useful bits (linear offset, direction vectors, parameters) and then format this nicely for TeX display. This is the only remaining TODO: add support for other vector input options, notably lists and ntuples. This will likely not happen in 2024. - Finished off `disp_eqns` and added `mat_disp_eqns` to do the same given `A` and `b` instead of a list of equations. - I've given up on PLU for now and have thus removed this section. I think that `rref`, `triangularize`, `echelon` and `lu_factor` do a good enough job, though it's still baffling to me that none of them do partial pivoting, nor do they avoid scaling rows unnecessarily.
Added test cases for v0.2.4
Added - rowscale and columnscale, intended as complements to rowop, rowswap, columnop and columnswap. It is odd that the third elementary row/column operation is not included by default. - setrowmx, setcolmx, and setdiagmx, intended as complements to setelmx. Used to overwrite rows, columns or diagonals.
Added tests for new functions in linearalgebra
There was a problem hiding this comment.
I've just had a quick pass through the functions you define. Mostly some nitpicking about names. (I find maxima not that consistent with its naming and some function names aren't very memorable, but maybe it's worth striving a bit for consistent naming.)
I will be authoring some linear algebra questions in the next month so will hopefully get to test this some more.
| * @param[list] ex A list of numbers | ||
| * @return[number] The greatest common divisor of all elements in ex | ||
| */ | ||
| lgcd(ex):= block([ex_gcd,ii], |
There was a problem hiding this comment.
You could probably simplify this by using lreduce (or similar) instead of a for loop: https://maxima.sourceforge.io/docs/manual/maxima_21.html#lreduce
| * @param[matrix] ta a mxn matrix | ||
| * @return[boolean] Are ex and ta row equivalent? | ||
| */ | ||
| row_equiv(ex,ta):= block( |
There was a problem hiding this comment.
Should this be called row_equivp, given that this is a predicate function?
Same for col_equiv and subspace_equiv.
There was a problem hiding this comment.
You're probably correct here. For some reason I had in my head that a predicate function would only take a single input, but a quick Google disproves that. I am wanting to do a (hopefully) final pass of function name changes before starting on question-writing in earnest, so I will fix these up then. My colleague had some other suggestions that I want to implement too.
| * @param[matrix] M a matrix | ||
| * @return[boolean] Is M a symmetric matrix? | ||
| */ | ||
| sym_p(M):= if squarep(M) then is(M = ev(transpose(M),simp)) else false; |
There was a problem hiding this comment.
The underscore in sym_p is inconsistent with the other predicate functions in its vicinity. Is symp already an existing keyword, or is there some other reason for this? Maybe symmetricp is a decent name too given that the others aren't abbreviated either?
There was a problem hiding this comment.
symmetricp already exists in the ctensor package but is forbidden, hence why I used the shortened name here as a maybe temporary option. I agree that symp is better than sym_p though, I'm not sure why I broke convention here. Perhaps I just read it is a misspelling of "simp" and didn't like it.
| * @param[matrix] M An mxn matrix | ||
| * @return[matrix] The same matrix with all entries below the diagonal set to 0 | ||
| */ | ||
| triu(M):= block([imax,jmax], |
There was a problem hiding this comment.
I get that diag is already taken as a keyword, but the naming inconsistency irks me a bit. Do you think it would be worthwhile renaming triu and tril to get_triu and get_tril?
There was a problem hiding this comment.
I'd rather keep triu and tril as they are and change get_diag to something else. I chose triu and tril to match the equivalent functions in NumPy because I am expecting to do a bit of Python CodeRunner alongside STACK. This is also why I chose column_stack as a function name. I wonder if something as simple as diagonal would work, it seems to be an unused function name.
| * @param[matrix] M An mxn matrix | ||
| * @return[boolean] Is this matrix upper triangular? | ||
| */ | ||
| triup(M):= block([tri,imax,jmax,ii,jj], |
There was a problem hiding this comment.
I noticed that you changed this from the earlier definition triup(M):= if matrixp(M) then is(M = triu(M)) else false; .
Just curious, is that for performance reasons or were there some other drawbacks?
There was a problem hiding this comment.
Whilst I was changing up triu etc I decided that it felt wrong to loop through the whole matrix to create triu(M) and then compare it element by element to M when I could just do it in one loop, so I changed it. I don't expect performance to be impacted massively (unless M is huge for some reason), but performance was indeed the motivation for the change.
Changed the names of some functions per suggestions from others; added two more vector parametric equation functions to test whether a given point is in a given affine subspace and whether a given vector is in a given vector subspace.
Adjusted test cases for recent commit to linearalgebra.mac
| * If v is a zero vector, returns false by default as the intended use case is determining whether a given DIRECTION vector is in a subspace. | ||
| * | ||
| * @param[matrix] v The vector of interest | ||
| * @param[list] dir_vecs A list of mx1 matrices that span the subspace. Does not need to be a basis. |
There was a problem hiding this comment.
I think when I wrote this myself, I chose to have a column matrix, where the columns spanned the subspace. I think with this choice it's more explicitly a "list of vectors", rather than an implicit convention on the columns of a matrix.
|
Thanks @LukeLongworth, that's a lot of work. Given it's in the "contrib" folder I'm happy to merge this into master to make it immediately available. I'm happy to do this. Are you ready? This is a lot of very useful work, thanks. |
|
I tested this on 4 different servers (our ecampus, our old sandbox, our new sandbox, and the Uni Trieste server). For 3 of them things run quite smoothly. The logs show: At the same time there are also unauthenticated requests to the Maxima service, which stop when cancelling the preview (closing the browser tab). For every 401 response from the Maxima service, there is a corresponding timeout message in the Moodle logs. It doesn't matter if the CAS connection timeout is set to 10 secs or 100 secs. This server uses goemaxima (like our other sandbox) and is hosted on google cloud platform, but has fewer resources: 'e2-small' (0.5-2 vCPU, 2 GB RAM) as opposed to 'e2-medium' (1-2 vCPU, 4 GB RAM). I don't think this is a Technically this is a resources/server configuration issue, and we can resolve it by e.g. giving our server more resources. But it does seem a bit scary that when importing a question from another server and trying to preview it, it could hang the whole server for a while (unless that's also a server setup issue -- seems a bit strange to me that a single request could hang the entire server). I wonder if it is worthwhile trying to modularize the library a bit to mitigate this risk, so that users don't have to include everything all the time? Is it possible to have recursive |
No. And in particular, it is impossible to conditionally include stuff. Having that inside the STACK compiler would be inconvenient, and I (as the one who built this first iteration of the include logic) don't believe it is worth it. Personally, I imagine that whatever is in the "contrib" will eventually end up in the core or become fragmented into small parts easier to include. I also imagine that at some point, someone will decide that they do not like this, and they build some wildly more complicated system based on the fact that the include just includes from an address. We might end up with someone setting up an external package management system which constructs a minimal include with all the required parts based on the URL, e.g.: |
|
I see, that makes a lot of sense. So with the vision of moving all or a lot of this into core, performance right now is not that relevant in the long run.
|
|
Yes, it ( Again, as long as we only work with URLs and have no special extra arguments for the inclusion outside that URL, one can simply tune the URL to get the version one wants. So that those who understand the risks and actually pay attention to these issues can just as well do something like this to get a very specific version of something: If someone wants to make something like |
|
Thanks for your testing of this @geoo89, I don't have the resources to do that myself and it sounds very useful. I don't think it's a big problem if I don't mind breaking this into smaller files, but would appreciate some guidance on best practice. My initial thoughts are:
I also agree with @aharjula's comments. I think that this project has stretched the limits of what @sangwinc I'm not finding myself making many more changes now (i.e. the occasional urge I'm getting to do more with it is more easily quashed) so you're welcome to add it to contrib if you like. As Georg mentioned, it might be more difficult to split up the file into smaller pieces whilst maintaining backwards compatibility. As a last comment/thought: if there are key functions that I think make sense as core maxima, where would those belong? Their own file alongside |
| if vectorp(ex) then return([list_matrix_entries(ex)]), | ||
| if is(op1="matrix") then return(args(transpose(ex))), | ||
| ex: args(ex), | ||
| if ev(every(lambda([ex2],constantp(ex2) or atom(ex2)),ex),simp) then return([ex]), |
There was a problem hiding this comment.
While working on some questions, I encountered some unexpected behavior:
I would expect make_list_of_lists([matrix([1],[2],[3]),matrix([2],[3],[4])]) to return [[1,2,3],[2,3,4]], but instead it returns
[[matrix([1],[2],[3]),matrix([2],[3],[4])]], i.e. wrapping the input into another layer of list brackets.
make_list_of_lists([c(1,2,3),c(2,3,4)]) also yields [[matrix([1],[2],[3]),matrix([2],[3],[4])]].
I observed that constantp(matrix([1],[2],[3])) evaluates to true, causing this early return. Some more investigation shows that constantp(matrix([1],[a],[3])) evaluates to false, and consequently make_list_of_lists([matrix([1],[2],[3]),matrix([2],[a],[4])]) and make_list_of_lists([c(1,2,3),c(2,a,4)]) return the expected output [[1,2,3],[2,a,4]].
This inconsistent behavior makes me believe that constantp/atom is maybe not the correct check to do, and maybe a similar type check for some type of iterable like in the first two lines is more appropriate?
Some more experimentation also shows that make_list_of_lists([1,a,g(x),4]) gives an error, because we don't do the early return here (neither constantp nor atom applies to g(x)) and the next line tries to call args on an atom.
As a more general question: What's the rationale for using lists of lists over lists of vectors (like in subspace_equivp and lin_indp) or a matrix (like lin_indp already supports)?
There was a problem hiding this comment.
Hi George,
I missed this message, sorry! I did push a fix to make_list_of_lists, as I had independently run into the same issue, and also identified that constantp wasn't doing what I expected. The new version just checks to see whether all elements in the list are scalars, and if so return a nested list. The rationale here is that if someone called make_list_of_lists([1,2,3]) in a linear algebra context I would personally expect the result to be [[1,2,3]] rather than [1,2,3] or [[1],[2],[3]]. I hope this change fixes the problems you highlighted!
To answer the more general question: Probably a coincidence, and I'm happy for extra functionality to be added here. Originally when I was writing these, I tried to account for every reasonable input type which meant that every function had a big list of if-thens at the beginning to put everything into a form that I could work with easily. Eventually I realised I was duplicating a lot of work and decided to move that "cleaning" process into a single function. My vague memory is that the lin_indp function supporting matrices is because I was expecting the use case of "does this matrix have full column rank" to be common enough that it would be annoying to use make_list_of_lists every time. Typing this explanation out has made me feel like I undercut my previous arguments somewhat!
Modify disp_eqns to support parameters
I didn't realise that vectors of constants were also considered constants, so one of the if-statements in make_list_of_lists wasn't behaving as intended. There may be a bug in point_in_affine_spacep; its 3rd test didn't work for me when running the whole bank of tests. I think that the problem is that I didn't declare local variables, so I have added that and I hope it fixes it.
%i is NOT the corresponding eigenvalue for matrix([0,-1],[1,0]) and c(1,%i), so the ta should be false. get_eigenvector returns an empty list rather than false when given an incorrect eigenvalue, so the ta should be [] instead of false.
|
Thanks for all your continued work on this @LukeLongworth. Rather than pull requests to master, should we set up a branch for this work. Then other people could contribute? |
|
@sangwinc that suits me fine! That could be a good place to experiment with splitting the file up into multiple smaller files too without interfering with the questions that myself and Georg are working on (I assume that we are both importing from my personal branch) |
|
Thanks @sangwinc ! I think a branch for contributions like this might be convenient. I've made some changes that are currently sitting on my fork (updated some functions and added a few). After some discussion with Luke I realized that the tests are meant to run with But we already observed that there is some functionality we are using a lot, like the |
|
@LukeLongworth, just looking through your library. It's a lot of work, and really impressive! Thank you for doing this. What would you like me to do next? |
- Fixed a bug in disp_eqns where negative pivots were missing their minus sign when the pivot was not in the first column. - Added tex1 into the format_as_coeff function as suggested by Georg Osang TODO: Unit tests need updating.
|
OK, we do need a branch for this work, otherwise were going to get conflicts between @LukeLongworth and @geoo89 and we don't need to do that! Luke, could you send a pull request here: https://github.com/maths/moodle-qtype_stack/tree/linear-algebra-beta @geoo89 I agree that |
As before, I don't intend for this pull request to be accepted yet. This is just me publishing the latest version for people to look at.
This is the latest version of #1148 which I will subsequently close. This version marks the point where I've stopped writing new functions (for now) and am now using this file to write my quizzes. I have included two files;
linearalgebra.macwhich is the file intended for eventual use, andlinearalgebra_no_test.macwhich is the same, but has nos_test_caselines, as these were causing issues when importing to questions withsimp:true. A separate issue, #1184, has been raised concerning this.The files contain:
Still to-do:
I think I've figured out my workflow now, and will keep everything in this branch, which I think means that I can update this pull request directly? I will figure that out when it comes to it.