I've started work on efficient sparse matrices over RDF/CDF, and here's my first round of questions/proposals. I commit to providing an implementation of these if accepted.
1) When multiplying sparse with dense, action.pyx converts sparse matrices to dense matrices. I'm not sure about other usecases, but for typical numerical usecases this is wrong -- if you have a sparse matrix it is typical to let it operate over several vectors at the same time (which would be expressed as a dense matrix). (In my usecases, I don't have enough memory (let alone time) for dense_matrix() to ever be called.) Proposal: Introduce _generic_matrix_times_matrix_, which can assume that both operands have the same parent *except* for sparsity. The default behaviour if sparsity is not the same is to ensure that both operands are dense, but this can be overridden. Disadvantage: Has to be overridden "on both sides", i.e. both for sparse and dense RDF/CDF matrices. Better ideas? 2) L = M.cholesky_decomposition() appears to return a writeable matrix. This causes some design problems for me. Permuting the matrix is essential to sparse Cholesky because a permutation can reduce the number of non-sparse entries introduced in L. Therefore L is typically (also with CHOLMOD, which I'll use) stored in a custom format. Assume one wants to perform L.solve(b). Then the options seem to be: a) Have cholesky_decomposition() return a custom class, Matrix_double_sparse_cholesky_factor, which can efficiently solve but is always immutable. __copy__() returns a writeable Matrix_double_sparse. b) Introduce M.solve_for_my_cholesky_factor(b). b) would lead to unnatural notation IMO, and if I want to avoid proxying-and-switch-backend-on-write, it seems that cholesky_decomposition() must return an immutable matrix by default. (In fact I'd wish for a convention that all matrix factorizations returned immutable matrices by default, since the factorizations are cached and one doesn't typically change them before using them.) 3) As mentioned earlier I'd like to implement explicitly diagonal and hermitian matrices (at least for RDF and CDF). Would this be OK?: sage: parent(hermitian_matrix(RDF, 3)) Full MatrixSpace of 3 by 3 Hermitian matrices over Real Double Field sage: parent(hermitian_matrix(CDF, 3)) Full MatrixSpace of 3 by 3 Hermitian matrices over Complex Double Field sage: parent(diagonal_matrix(RDF, 3, enforce=True)) Full MatrixSpace of 3 by 3 diagonal matrices over Real Double Field It could kind of work to simply add flags about sparsity patterns and keep them in the sparse matrixspace, but ideally I'd like them to raise exceptions when assigning outside the lower triangle/diagonal, and that seems to call for another parent? Also these spaces would only allow square matrices I think. Dag Sverre -- To post to this group, send an email to sage-devel@googlegroups.com To unsubscribe from this group, send an email to sage-devel+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-devel URL: http://www.sagemath.org
