Lattice points at the intersection of 3 hyperplanes

[mrKirosana]

I'm wondering whether polymake can help me with the following problem.

I wish to enumerate all lattice points contained in the intersection of three hyperplanes in the non-negative half-space. That is to say, I have three linear equations of d variables and I would like to find all non-negative integer solutions.

I guess one should be able to consider the above as a polyhedron where the faces are the three hyperplanes and the ones arising from restricting to non-negative solutions only, but would the lattice points at the intersection be a subset of the vertices in that case? If so, would there be any way to find the subset directly, or is the only option to find all vertices and then filter?

[joswig]

If you replace what obviously needs to be replaced in

Code: Select all

$p = new Polytope(INEQUALITIES=>{non-negativity constraints}, EQUATIONS=>{three hyperplanes}); print $p->LATTICE_POINTS;

then this should do the trick.

Several caveats: (1) depending on your hyperplanes the result may be unbounded. (2) Deciding whether a polyhedron contains one lattice point or not is an NP-hard problem. Therefore this will not work in higher dimensions. What exactly "higher" means cannot be said easily.

If you want to learn more, please check out the following:

Assarf et al.: Computing convex hulls and counting integer points with polymake, Mathematical Programming Computation,
March 2017, Volume 9, Issue 1, pp 1-38
https://link.springer.com/article/10.10 ... 016-0104-z

[opfer]

I wish to enumerate all lattice points contained in the intersection of three hyperplanes in the non-negative half-space. That is to say, I have three linear equations of d variables and I would like to find all non-negative integer solutions.

I would be happy if you could share your equations with us. I am currently working on an additional algorithm for such problems and could test it on your data.

Best regards,
Thomas

[mrKirosana]

I wish to enumerate all lattice points contained in the intersection of three hyperplanes in the non-negative half-space. That is to say, I have three linear equations of d variables and I would like to find all non-negative integer solutions.

I would be happy if you could share your equations with us. I am currently working on an additional algorithm for such problems and could test it on your data.

Best regards,
Thomas

The equations are in fact of rather specific form. Namely, you consider a range of positive integers from \( i_{min} \) to \( i_{max} \), and positive integer constants \( N \), \( D \) and \( A \). Then you may write

\( m.x=b \),

where the rows of m are as follows. First row is \( i_{min}^2, (i_{min}+1)^2,...,i_{max}^2 \), second row is \( i_{min}, i_{min}+1,...,i_{max} \), and the final row is all 1s. Column vector \( x \) contains the variables and column vector \( b \) reads \( A \), \( D \), \( N \). I am looking for all non-negative integer solutions of \( m.x=b \).

For a small concrete example, consider \( i_{min}=1 \), \( i_{max}=3 \), and \( b=(70,32,15)^T \). In this particular case there is only one solution and it is non-negative, namely \( (3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2)^T \).

Due to the peculiar form of matrix \( m \), one approach to enumerating the solutions is to find the integer partitions of \( A \) into exactly \( N \) elements using only \( i_{min}^2, (i_{min}+1)^2,...,i_{max}^2 \). Once you have the partititions, you can then take the element-wise square root of each and select the ones that add up to \( D \). This is not too bad, I have used it to solve cases where for example \( b=(1284,168,30)^T \), \( i_{min}=2 \), \( i_{max}=17 \). On the other hand, in this bigger example only 9 % of the found partitions add up to \( D \) so while this works it is wasteful, and so I have been looking for a more direct method.

[opfer]

Thank you for the information. I ran my code on both examples.

For the first example, I think you provided the wrong solution. The example is three-dimensional and your solution is 15-dimensional. I think the correct solution is \( (0,13,2)^T \). This is found instantaneously (for obvious reasons).

I also ran your second example for two hours. It found a lot of lattice points, but it did not terminate. How long does your approach take? I guess, it is much faster as it is a specialized algorithm for your problem, not a general algorithm for general polytopes.

Do you have good examples for A, D and N in lower dimensions?

Best regards,
Thomas

[mrKirosana]

Thank you for the information. I ran my code on both examples.

For the first example, I think you provided the wrong solution. The example is three-dimensional and your solution is 15-dimensional. I think the correct solution is \( (0,13,2)^T \). This is found instantaneously (for obvious reasons).

I also ran your second example for two hours. It found a lot of lattice points, but it did not terminate. How long does your approach take? I guess, it is much faster as it is a specialized algorithm for your problem, not a general algorithm for general polytopes.

Do you have good examples for A, D and N in lower dimensions?

Best regards,
Thomas

Right, for the first example \( (0,13,2)^T \) is indeed the expected output. These numbers can basically be interpreted as multiplicities for the elements appearing in the matrix \( m \), so what I said earlier is the result of "take 0 1s, 13 2s and 2 3s".

Your news about the second example are a bit disheartening, as my code takes around 15 seconds for that. It's basically a one-liner with Wolfram Mathematica, based on the integer partitioning approach. The problem with that, besides being wasteful, is that it is very memory-intensive and larger examples will simply crash my computer. Perhaps I should look for efficient integer partitioning algorithms implemented in some low level language..

A straightforward way to generate more \( m \) and \( b \) to try is to consider degree sequences of random graphs, and set \( N \) to be the number of nodes, \( D \) the total degree, \( A \) the sum of squared degrees and set \( i_{min} \), \( i_{max} \) to minimum and maximum degree, respectively. This of course guarantees that there will always be at least one solution. I generated a random tree of 30 nodes, and got \( b=(148,58,30)^T \) with \( i_{min}=1 \), \( i_{max}=6 \). This gave me 14 solutions.

[opfer]

Hi,

I re-ran my code with your second example and it terminated after 3-4 hours. It found 287301 lattice points. Is this the same number that you get?

My approach is slow, because it does a lot of operations on rational coordinates. I am not very surprised that a specialized algorithm is much faster. At least, my algorithm does not consume too much memory.

You might try the other algorithms that are implemented in polymake, but I would not expect too much. If you have runtimes of 15 seconds for your code, maybe this is fast enough? Maybe, you can find out why it uses so much memory.

I will also play around with your new example.

Best regards,
Thomas

[opfer]

I generated a random tree of 30 nodes, and got \( b=(148,58,30)^T \) with \( i_{min}=1 \), \( i_{max}=6 \). This gave me 14 solutions.

For this example, I also get 14 solutions (instantaneously).

Best regards,
Thomas

[mrKirosana]

I re-ran my code with your second example and it terminated after 3-4 hours. It found 287301 lattice points. Is this the same number that you get?

Yes.

You might try the other algorithms that are implemented in polymake, but I would not expect too much. If you have runtimes of 15 seconds for your code, maybe this is fast enough? Maybe, you can find out why it uses so much memory.

Speed is not a big issue here, but the memory consumption is. The key part of my code (Wolfram Language) is simply

Code: Select all

IntegerPartitions[restrictionA, {restrictionN},Range[dmin, dmax]^2]]

so I can't really tweak this at all, internally it is some integer partitioning algorithm implemented in C. For the bigger example of 287301 the intermediate expression returned by the above line is around 2.5 GB while the final result, i.e. the lattice points, only take around 70 MB.

Basically I was hoping to experiment with the parameter \( N \) going up to 100, but it's starting to look pretty bad. Of course, \( N \) alone is not the only thing that matters. If you take as a starting point the degree sequence of a regular graph, i.e. all degrees are some constant, then even \( N=600 \) evaluates at once. The computational complexity depends heavily on the combination of all parameters it seems.

EDIT

You can play around with limited size examples at Wolfram cloud sandbox, just copy the code below. To get started, you need to make a new input cell by clicking the circled + and choosing "Wolfram Language input". Evaluation is done with shift-enter. Here the output would be the number of solutions.

Code: Select all

{restrictionN,restrictionD,restrictionA,dmin,dmax}={30,58,146,1,5}; solutions = Sqrt@Select[IntegerPartitions[restrictionA, {restrictionN}, Range[dmin, dmax]^2], Total[Sqrt[#]] == restrictionD &]; Length@solutions

[opfer]

I re-ran my code with your second example and it terminated after 3-4 hours. It found 287301 lattice points. Is this the same number that you get?

Yes.

Good. That indicates that my algorithm works correctly

You might try the other algorithms that are implemented in polymake, but I would not expect too much. If you have runtimes of 15 seconds for your code, maybe this is fast enough? Maybe, you can find out why it uses so much memory.

Speed is not a big issue here, but the memory consumption is. The key part of my code (Wolfram Language) is simply

Code: Select all

IntegerPartitions[restrictionA, {restrictionN},Range[dmin, dmax]^2]]

so I can't really tweak this at all, internally it is some integer partitioning algorithm implemented in C. For the bigger example of 287301 the intermediate expression returned by the above line is around 2.5 GB while the final result, i.e. the lattice points, only take around 70 MB.

Basically I was hoping to experiment with the parameter \( N \) going up to 100, but it's starting to look pretty bad. Of course, \( N \) alone is not the only thing that matters. If you take as a starting point the degree sequence of a regular graph, i.e. all degrees are some constant, then even \( N=600 \) evaluates at once. The computational complexity depends heavily on the combination of all parameters it seems.

I guess Mathematica uses some dynamic programming approach storing lots of intermediate results to solve this problem so quickly. At least this would explain why it consumes so much memory. I have no idea whether there is an algorithm that is fast and has a small memory footprint for your problem.

Do you have a machine with more memory available?

Best regards,
Thomas