4
$\begingroup$

The kernel to the left of the matrix

Nm = {{
    -1, 1, -1, 1, 1, 0},
   {2, -2, 0, 0, -1, 1 },
   {0, 0, -1, 1, 1, 0 },
   {0, 0, 1, -1, 0, -1},
   {0, 0, 0, 0, -1, 1}};

may be generated by {0, 0, 1, 1, 1},{2, 1, 0, 2, 1}, as easily checked

{0, 0, 1, 1, 1} . Nm
{2, 1, 0, 2, 1} . Nm

I need a script which finds a positive (nonnegative) basis for arbitrary integer valued matrix. For now, I am working with

tx = Table[x[i], {i, 5}];Solve[Thread[tx . Nm == 0], tx]

and then searching with pencil for an answer. I tried things like adding positivity conditions

Solve[Join[Thread[tx . Nm == 0], Thread[tx >= 0]], tx]

or

FindInstance[Join[Thread[tx . Nm == 0], Thread[tx > 0]], tx, 2]

but they don't seem to advance me, so I'm asking for help.

Here's a second example, where all the methods offered are in difficulty

ST = {{-1, 0}, {0, 1}, {1, -1}, {-2, 1}, {1, -1}}
nu=NullSpace[ST // Transpose]

The null-space contains two positive vectors out of three ({{1, 1, 0, 0, 1}, {-2, -1, 0, 1, 0}, {1, 1, 1, 0, 0}})

The

HermiteDecomposition[nu][[2]]
(*{{1, 0, 0, -1, -1}, {0, 1, 0, 1, 2}, {0, 0, 1, 0, -1}}*)

discovers a third vector. The script I offered has 5 arbitrary constants instead of three, and after playing around with substituting values and reducing size yields another answer {0,1,2,1,0}. The impressive algorithm of @Daniel Lichtblau does not seem to strive for minimality. So, the question of offering a script for this seemingly elementary problem remains open :)

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ "they don't seem to advance me" - what's the problem exactly? Do you get solutions? Try xx = Array[x, 5]; FindInstance[{xx . Nm == 0, xx >= 0}, xx, Reals, 3]. Are those solutions good for your purposes? You could also skip the explicit positivity condition and replace it with the PositiveReals domain specification. $\endgroup$
    – MarcoB
    Commented Dec 7, 2023 at 11:47
  • $\begingroup$ FindInstance[{xx . Nm == 0, xx >= 0}, xx, Integers, 3] almost does it. I only need to divide the answers {x[1] -> 226, x[2] -> 113, x[3] -> 0, x[4] -> 226, x[5] -> 113}, {x[1] -> 0, x[2] -> 0, x[3] -> 12, x[4] -> 12, x[5] -> 12} by 113, and 12. But I still would like to add a condition that the GCD of each solution is 1, and I'm not sure how to code that $\endgroup$
    – florin
    Commented Dec 7, 2023 at 12:08
  • $\begingroup$ I can't even use GCD, since it does not accept lists as arguments, but only "enumerations" of arguments, and I don't know how to convert a list into an enumeration: xx = Array[x, 5] x = xx /. FindInstance[{xx . Nm == 0, xx >= 0}, xx, Integers, 3] GCD[Flatten[x[[2]]]] does not work $\endgroup$
    – florin
    Commented Dec 7, 2023 at 12:25
  • $\begingroup$ FindInstance does not guarantee independence, as illustrated by the second example I added $\endgroup$
    – florin
    Commented Dec 7, 2023 at 12:57
  • $\begingroup$ If you have a list m = {3, 6, 9, 12}; then you can Apply the GCD function to that list as follows: GCD @@ m. That's a general solution to the problem of using functions that accept sequences of arguments on lists. $\endgroup$
    – MarcoB
    Commented Dec 8, 2023 at 3:22

5 Answers 5

Reset to default
4
$\begingroup$

At least for your example, NullSpace does precisely what you want,

NullSpace[Fm]

(* {{1, 0, 1, 1, 0, 1}, {0, 0, 0, 1, 1, 0}, {1, 1, 0, 0, 0, 0}} *)

Actually I just realized that what you want is the co-kernel instead, the following seems to do the trick:

NullSpace[Fm // Transpose]
Array[c, Length[%]] . %
FindInstance[Flatten@{Thread[Array[c, Length[%%] {1, 1}] . %% >= 0], Det[Array[c, Length[%%] {1, 1}]] != 0}, Flatten@Array[c, Length[%%] {1, 1}], Integers]
ArrayReshape[Last /@ %[[1]], Length[%%%] {1, 1}] . %%%

(* {{0, 1, 1, 1, 0, 1}, {1, 0, 1, 0, 0, 0}, {0, 0, 0, 0, 1, 1}} *)
Improve this answer
$\endgroup$
5
  • $\begingroup$ Sorry, on my computer (with 13.3) the code below yields mistakes, and I do not understand what it is supposed to do Nm = {{ -1, 1, -1, 1, 1, 0}, {2, -2, 0, 0, -1, 1 }, {0, 0, -1, 1, 1, 0 }, {0, 0, 1, -1, 0, -1}, {0, 0, 0, 0, -1, 1}}; NullSpace[Nm // Transpose] Array[c, 3] . % FindInstance[ Flatten@{Thread[Array[c, {3, 3}] . %% >= 0], Det[Array[c, {3, 3}]] != 0}, Flatten@Array[c, {3, 3}], Integers] ArrayReshape[Last /@ %[[1]], {3, 3}] $\endgroup$
    – florin
    Commented Dec 8, 2023 at 6:50
  • $\begingroup$ It seems FindInstance cannot be part of a solution, since it does not produce algebraically independent vectors $\endgroup$
    – florin
    Commented Dec 8, 2023 at 6:52
  • $\begingroup$ NullSpace[Fm] produces vectors with mixed signs $\endgroup$
    – florin
    Commented Dec 8, 2023 at 6:53
  • $\begingroup$ @florin You need to replace 3 by the dimension of the co-kernel, I updated the code to make this automatic. $\endgroup$
    – AccidentalFourierTransform
    Commented Dec 8, 2023 at 16:08
  • $\begingroup$ Thanks, it's a quite nice direct solution, but not easy to decipher. Maybe writing like nul = NullSpace[Fm // Transpose]; dim = nul // Length; FindInstance[ Flatten@{Thread[Array[c, dim {1, 1}] . nul >= 0], Det[Array[c, dim {1, 1}]] != 0}, Flatten@Array[c, dim {1, 1}], Integers]; lincom = ArrayReshape[Last /@ %[[1]], dim {1, 1}]; lincom . nul would have made it easier. And I must confess I could not decipher the line before the last $\endgroup$
    – florin
    Commented Dec 11, 2023 at 10:36
3
$\begingroup$ 
Nm = {{-1,  1, -1,  1,  1,  0},
      { 2, -2,  0,  0, -1,  1},
      { 0,  0, -1,  1,  1,  0},
      { 0,  0,  1, -1,  0, -1},
      { 0,  0,  0,  0, -1,  1}};

X = Array[x, Length[Nm]];
sol = SolveValues[Thread[X . Nm == 0], X, NonNegativeIntegers]

$$ \{\{ 2 c_1\text{ if }(c_1|c_2)\in \mathbb{Z}\land c_1\geq 0\land c_2\geq 0,\\ c_1\text{ if }(c_1|c_2)\in \mathbb{Z}\land c_1\geq 0\land c_2\geq 0,\\ c_2\text{ if }(c_1|c_2)\in \mathbb{Z}\land c_1\geq 0\land c_2\geq 0,\\ 2 c_1+c_2\text{ if }(c_1|c_2)\in \mathbb{Z}\land c_1\geq 0\land c_2\geq 0,\\ c_1+c_2\text{ if }(c_1|c_2)\in \mathbb{Z}\land c_1\geq 0\land c_2\geq 0 \}\} $$

By picking specific values for the coefficients $(c_1,c_2)\in\mathbb{N}_0$ we get non-negative integer vectors:

sol /. {C[1] -> 1, C[2] -> 0}
(*    {{2, 1, 0, 2, 1}}    *)

sol /. {C[1] -> 0, C[2] -> 1}
(*    {{0, 0, 1, 1, 1}}    *)
Improve this answer
$\endgroup$
2
  • $\begingroup$ im worried that you never impose linear independence, so it is not clear that this will always yield a basis of the co-kernel. You might have to try lots of different choices for C until one gives you a basis, right? $\endgroup$
    – AccidentalFourierTransform
    Commented Dec 13, 2023 at 14:17
  • $\begingroup$ @AccidentalFourierTransform yes, in general that's right. $\endgroup$
    – Roman
    Commented Dec 13, 2023 at 18:40
1
$\begingroup$

A method that is not guaranteed but might be a reasonable heuristic is to take the Hermite form of the null vectors. This guarantees that the first nonzero elements are all in different positions, and all positive. So at worst you might need to take nonnegative combinations of these to get all-nonnegative vectors. I'll say a bit more about this below.

fm = {{-1, 1, 1, 0, 0, 0}, {-1, 1, 0, 0, 0, 1}, {1, -1, -1, 0, 0, 
    0}, {0, 0, 1, -1, 1, 0}, {0, 0, 0, -1, 1, 1}, {0, 0, 0, 
    1, -1, -1}};
nullsfm = HermiteDecomposition[NullSpace[Transpose[fm]]][[2]]

(* Out[870]= {{1, 0, 1, 0, 0, 0}, {0, 1, 1, 1, 0, 1}, {0, 0, 0, 0, 1, 1}} *)

Another method is guaranteed except it might run forever. It uses integer programming in a way similar in spirit to the ones already shown. we start by finding a null space basis. We create a transformation matrix and enforce that it be unimodular (det = +-1) and that the new vectors be nonnegative and have at least one nonzero element.

fm = {{-1, 1, 1, 0, 0, 0}, {-1, 1, 0, 0, 0, 1}, {1, -1, -1, 0, 0, 
    0}, {0, 0, 1, -1, 1, 0}, {0, 0, 0, -1, 1, 1}, {0, 0, 0, 
    1, -1, -1}};
nullsfm = NullSpace[Transpose[fm]]
n = Length[nullsfm];

(* Out[872]= {{-1, 1, 0, 1, 0, 1}, {1, -1, 0, -1, 1, 0}, {1, 0, 1, 0, 0, 0}} *)

aa = Array[a, {n, n}];
faa = Flatten[aa];
prod = aa . nullsfm;
c1 = Thread[Flatten[prod] >= 0];
c2 = Thread[Map[Total, prod] >= 1];
constraints = Join[c1, c2, {Det[aa]^2 == 1}];
{min, vals} = Minimize[{1, constraints}, faa, Integers]

(* Out[886]= {1, {a[1, 1] -> 1, a[1, 2] -> 0, a[1, 3] -> 1, a[2, 1] -> 0,
   a[2, 2] -> 0, a[2, 3] -> 1, a[3, 1] -> 1, a[3, 2] -> 1, 
  a[3, 3] -> 0}} *)

Get the new nulls.

posnulls = prod /. vals

(* Out[902]= {{0, 1, 1, 1, 0, 1}, {1, 0, 1, 0, 0, 0}, {0, 0, 0, 0, 1, 1}} *)

Check that they are in fact nulls.

posnulls . fm

(* Out[903]= {{0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}} *)

This method is likely to be persnickety in that some minor variants seem to hang.

I prefer the first method though. It seems like it can be completed to an algorithm without resorting to integer programming, so the time complexity would then be guaranteed to be reasonable. At worst it would require integer linear programming to determine nonnegative combinations of each other null vector to add to a given one. Which might be equivalent to the solution you show, I'm not sure.

Improve this answer
$\endgroup$
6
  • $\begingroup$ This problem is very importaht for the determination of the conservation laws of interaction network ODE's, but the actual computational implementation falls outside the scope of interest of the literature (too simple for some, too hard for others?). I feel a relief now that stackexchange has filled this gap $\endgroup$
    – florin
    Commented Dec 12, 2023 at 9:03
  • $\begingroup$ This is to find basis vectors for a compatibility class in the positive orthant? You might be correct about the "too simple...too hard" assessment. $\endgroup$
    – Daniel Lichtblau
    Commented Dec 12, 2023 at 16:21
  • $\begingroup$ Yes, it's a fundamental step both for finding basis vectors for a compatibility class in the positive orthant , and for finding positive basis vectors in the kernel to the left (these yield conservation laws) $\endgroup$
    – florin
    Commented Dec 12, 2023 at 18:02
  • $\begingroup$ is it important that the matrix is unimodular (rather than being just invertible over the rationals)? I guess you want a basis that generates a lattice that is in a sense minimal, i.e., without (integer) sublattices? this seems useful but perhaps not important for the OP, and dropping this condition will likely speed up the search by a lot, no? $\endgroup$
    – AccidentalFourierTransform
    Commented Dec 13, 2023 at 19:10
  • $\begingroup$ @AccidentalFourierTransform You cannot entirely get away from the extra condition insofar as the conversion matrix must have full rank. $\endgroup$
    – Daniel Lichtblau
    Commented Dec 13, 2023 at 21:36
1
$\begingroup$

The script I needed is a combination of the two answers offered by @Roman and @AccidentalFourierTransform:

mat = {{-1, 1, 1, 0, 0, 0}, {-1, 1, 0, 0, 0, 1}, {1, -1, -1, 0, 0, 
    0}, {0, 0, 1, -1, 1, 0}, {0, 0, 0, -1, 1, 1}, {0, 0, 0, 
    1, -1, -1}};
con[mat_] := Module[{X, sol, dim, cv}, X = Array[x, Length[mat]];
  sol = SolveValues[Thread[X . mat == 0], X, NonNegativeIntegers];
  dim = NullSpace[mat // Transpose] // Length;
  cv = Table[C[i], {i, dim}];
  Flatten /@ 
   Table[sol /. Thread[cv -> IdentityMatrix[dim][[i]]], {i, dim}]]
con[mat]
Improve this answer
$\endgroup$
1
$\begingroup$

My first answer was not really adequate. Here is another way, related and also heuristic, but I think more likely to succeed. We again get the Hermite form (HNF). We use a bit if integer programming to get an integer combination from these that is nonnegative and furthermore strictly positive in all columns that do not have pivots in the HNF. For code purposes I take a shortcut here and also omit finding the transformation matrix, and furthermore assume there are no zero columns. All these simplifications are easy to redress.

We form this nonnegative vector with the added stipulation that the first vector in the HNF is multiplied by 1. This will allow us to enforce the condition that the conversion matrix be invertible over the integers. The heuristic part is that I have no proof this can always be done. (So maybe it's obvious, maybe it's wrong...) Were it to fail one could of course use a different basis vector from the HNF, but again, I don't know if one can always get a positive vector in this way.

With this vector we can form a full basis that is nonnegative. I'll illustrate with an example.

First I start with a set of random linearly independent nonnegative integer vectors.

SeedRandom[111123];
nvecs = 7;
dim = 12;
randPosVecs = RandomInteger[{0, 10}, {nvecs, dim}];
MatrixRank[posvecs] == nvecs

(* Out[1096]= True *)

I use a utility to create a "random" unimodular matrix by multiplying an upper triangular matrix with ones on the main diagonal with a lower triangular matrix similarly constructed.

randUnimodularMatrix[n_] := Module[
  {top, bot, bands},
  top = IdentityMatrix[n] + 
    SparseArray[
     Table[Band[{1, j}] -> RandomInteger[{-10, 10}, n - j + 1], {j, 2,
        n}], {n, n}];
  bot = IdentityMatrix[n] + 
    SparseArray[
     Table[Band[{j, 1}] -> RandomInteger[{-10, 10}, n - j + 1], {j, 2,
        n}], {n, n}];
  top . bot]

We'll create one and check unimodularity.

randuni = randUnimodularMatrix[nvecs];
Det[randuni]

(* Out[1098]= 1 *)

Now we apply this to our nonnegative basis.

vecs = randuni . randPosVecs

(* Out[1110]= {{-28, -207, -276, -359, -53, -192, -446, -374, -252, \
-377, -243, -178}, {274, -54, 242, 297, 253, 285, 321, 495, 96, 415, 
  795, 276}, {29, -293, -302, -461, -265, -380, -621, -754, -102, \
-587, -426, -516}, {42, -450, -464, -679, -302, -493, -913, -905, \
-248, -741, -433, -617}, {279, 114, 193, -374, -55, 113, -334, -648, 
  192, -361, -749, -197}, {-35, -89, -98, -32, -30, -88, -73, -18, \
-69, -44, 78, -54}, {3, 84, 81, 72, 36, 80, 111, 90, 51, 86, -7, 79}} *)

So we have a mix of positive and negative values. Furthermore the HNF is not all nonnegstive. We will use it to create a new basis though.

newvecs = HermiteDecomposition[vecs][[2]]

(* Out[1111]= {{1, 0, 0, 0, 0, 0, 6681, 10707, -9738, 8753, 14524, 
  8735}, {0, 1, 0, 0, 0, 0, 19558, 31395, -28519, 25667, 42568, 
  25591}, {0, 0, 1, 0, 0, 0, 13724, 22079, -20023, 18053, 29923, 
  17976}, {0, 0, 0, 1, 0, 0, 13508, 21689, -19697, 17732, 29407, 
  17676}, {0, 0, 0, 0, 1, 0, 20948, 33657, -30554, 27517, 45627, 
  27422}, {0, 0, 0, 0, 0, 1, 15135, 24311, -22073, 19876, 32959, 
  19810}, {0, 0, 0, 0, 0, 0, 22488, 36125, -32798, 29535, 48975, 
  29435}} *)

Here is the ILP step where we create a vector that is noonnegative, with strictly positive values in all non-pivot columns.

tvars = Array[t, nvecs];
t[1] = 1;
posvec = tvars . newvecs;
ptot = Total[posvec];
{min, vals} = 
  Minimize[{ptot, 
    Join[Thread[posvec >= 0], 
     Thread[posvec[[Length[newvecs] + 1 ;;]] >= 1]]}, Rest@tvars, 
   Integers];
tvals = tvars /. vals;
posvec = tvals . newvecs

(* Out[1125]= {1, 3, 4, 9, 11, 0, 3, 11, 4, 6, 5, 3} *)

Now we create a full basis. I add appropriate multiples of posvec to the second through last vector in newvecs, which we know is a basis set. I use posvec itself for the first basis vector. It is straightforward to show from the construction that this is in fact a basis set, that is, comes from an invertible (over the integers) transformation of the original basis.

nnvecs = Table[
   vec = newvecs[[j]];
   qvec = 
    Table[If[posvec[[k]] == 0, 0, Floor[vec[[k]]/posvec[[k]]]], {k, 
      Length[vec]}];
   qmin = Min[qvec, 0];
   vec - qmin*posvec
   , {j, Length[newvecs]}
   ];
nnvecs[[1]] = posvec;

First check that these are all nonnegative.

Map[MinMax, nnvecs]

(* Out[1081]= {{0, 9}, {1, 13863}, {0, 10558}, {0, 111614}, {0, 36362}, {0, 40519}, {0, 185970}} *)

Next show that this basis comes from an invertible transformation of newvecs (so all integers and determinant is +-1).

cmat = LinearSolve[Transpose@newvecs, Transpose@nnvecs];
Union[Flatten[Map[IntegerQ, cmat, {2}]]]
Det[cmat]

(* Out[1087]= {True}

Out[1088]= 1 *)

Though this is redundant, we also show that it is a basis for the original vecs by noting that they have the same HNF.

HermiteDecomposition[nnvecs][[2]] === HermiteDecomposition[newvecs][[2]]

(* Out[1090]= True *)

So the open question is whether this can be made into a guaranteed algorithm.

Also on the unfortunate side is that the result, while a basis, is not going to deliver what I believe are called extremal or elementary fluxes (in reference to the underlying problem, which I believe is to do compoutations with a mass-action chemical reaction network).

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.