John Cremona on Thu, 21 Nov 2024 16:33:09 +0100 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: integral homology |
Dear John,
On 21/11/2024 14:45, John Cremona wrote:
Would you be interested in trying my experimental code? You would have to convert your input matrices into the Pari representation of sparse matrices (zMs), namely a t_VEC of columns zCs, each being a t_VEC with two components, a t_VECSMALL of indices and a t_VECSMALL of coefficients. For instance:On Thu, 21 Nov 2024, 13:00 Aurel Page, <aurel.page@normalesup.org> wrote:On 21/11/2024 12:40, John Cremona wrote:
This is fine for the algorithm I mentioned, but only the invertible coefficients will be used as pivots in the sparse elimination, so the algorithm will be more useful if there are many +-1 among the nonzero coefficients.Thanks, Aurel and Watson. As either of you might have guessed, my context is integral homology of Bianchi orbifolds, hyperbolic 3 space modulo GL2 or SL2 if the ring of integers of an imaginary quadratic field. The matrix is very sparse but entries might be bigger than 1, eg the boundary if a triangle could be 3 times the same edge.
I think the vast majority will be +-1.
By the way, the size of the matrices currently being processed is about 25000 x 15000
M = [5,0,2; 0,0,-1; 4,7,0]
is represented as
[[Vecsmall([1,3]),Vecsmall([5,4])], [Vecsmall([3]),Vecsmall([7])], [Vecsmall([1,2]),Vecsmall([2,-1])]]
Best,
Aurel
On Thu, 21 Nov 2024, 10:22 Aurel Page, <aurel.page@normalesup.org> wrote:Dear John,
There are several questions that influence the answer:
1/ do you only want the structure of the homology or also a basis?
2/ are your matrices sparse?
If you want a basis, the way I have implemented it is a variant of what you proposed:
Z = matker(D1); \\Q-kernel
Z = matrixqz(Z,-2); \\Z-kernel
if(n>0 && matsize(Z)[1]!=0,
D2 = matconcat([D2,vectorv(n)]); \\trick for matrices with 0 columns :-(
);
B = matsolve(Z,D2);
matsnf(B,4);
However, if you only want the structure, you don't need to compute the kernel. The structure of the torsion in the homology will be directly given by the SNF of the second differential, while the rank of the homology is simply the dimension of the middle space minus the sum of the ranks of the two differentials, so you only need one SNF and one rank computation, no integer kernel.
About sparseness: I have been playing around with an algorithm for computing homology of a complex with sparse matrices with most coefficients +-1 (the core algorithm is in the branch aurel-amt, but there is no documentation and it is still work in progress). It has been working well on my examples. For instance, the homology of an arithmetic hyperbolic 5-orbifold:
dimensions: [0, 15151, 251927, 1079420, 1887244, 1462440, 417840, 0]
number of nonzero coefs: [0, 503854, 3509856, 8426440, 8356800, 2924880, 0]
The structure of the entire homology over Z was computed in 5 minutes with 6.4GB of memory (single core on my laptop).
Best,
Aurel
On 18/11/2024 10:02, John Cremona wrote:
Well, I no longer need an answer to my question since it turned out that in my specific application, the matrix A10 was so simply structured that I could simply replace all of steps 1-3 by a new preliminary step which just deleted certain rows from A21 without having to do the HNF or matrix inversion steps.
It would still be nice to know if there is a better way to do the general case, For example, in the HNF computation of H and U from A10, I presume that a series of column operations are applied to A10, with U recording those operations; and then it would be easy to compute U^{-1} at the same time step by step, I think. In some of the examples I was having difficulty with where the dimensions n0, n1, n2 are 36,14219, and 24599, I was running out of memory even after setting PARISIZEMAX to 50000000000, i.e. 50G, and this occurred while in _ZM_inv_worker() so while inverting U.
John
On Fri, 15 Nov 2024 at 15:03, John Cremona <john.cremona@gmail.com> wrote:
I am computing some integral homology in a C++ program using libpari, and while what I have works fine it is possible that there is a better (simpler or faster) way to do this, possibly even a built-in function which I have not found which does what I want more simply.
I have Z-linear maps maps Z^n2 -> Z^n1 -> Z^n0 whose composite is 0 with integer matrices A10 (size n0 x n1) and A21 (size n1 x n2), so acting on column vectors on the right, and A10*A21=0. All I want is the abelian group structure of ker(A10)/im(A21). A rough sketch of what I do is this:
1. Find the HNF of A10, say H = A10*U with U in GL(n1,Z) representing a change of basis for the module in the middle.2. Compute M = U^{-1}*A21, the new matrix of the second map.3. Drop the last r rows of M, where r is the rank of H (also of A10).4. Return the SNF of M.
Is there a better way?
John