Simultaneous linear solving
Simultaneous linear solving
Set up: Create an example with known common solution
Set up: Create an example with known common solution
SeedRandom[1234];mat1=RandomInteger[{-5,5},{3,3}]mat2=RandomInteger[{-5,5},{3,3}]solc=RandomInteger[{-5,5},3]vec1=mat1.solcvec2=mat2.solc
Out[]=
{{-5,1,4},{1,5,-5},{2,-5,-5}}
Out[]=
{{3,-1,-1},{3,0,4},{2,-3,3}}
Out[]=
{-1,0,3}
Out[]=
{17,-16,-17}
Out[]=
{-6,9,7}
Part 1: Solve first matrix/vector
Part 1: Solve first matrix/vector
newmat1=Join[-{vec1},Transpose@mat1];augmat1=Join[newmat1,IdentityMatrix[Length[newmat1]],2]
Out[]=
{{-17,16,17,1,0,0,0},{-5,1,2,0,1,0,0},{1,5,-5,0,0,1,0},{4,-5,-5,0,0,0,1}}
Row reduce. Row with 0s in first three slots and a 1 in the fourth has solution values after that 1.
In[]:=
redmat1=RowReduce[augmat1]soln=FirstCase[redmat1,{Repeated[0,{Length[vec1]}],1,sol__}:>{sol}]
Out[]=
1,0,0,0,-,-,-,0,1,0,0,-,,-,0,0,1,0,-,-,-,{0,0,0,1,-1,0,3}
10
37
1
37
3
37
3
37
17
185
23
185
5
37
21
185
26
185
Out[]=
{-1,0,3}
Part 2: Find common solution for both matrix/vector pairs
Part 2: Find common solution for both matrix/vector pairs
In[]:=
bigvec=Join[vec1,vec2];bigmat0=Apply[Join[##,2]&,Map[Transpose,{mat1,mat2}]];bigmat=Join[-{bigvec},bigmat0];augmat=Join[bigmat,IdentityMatrix[Length[bigmat]],2]
Out[]=
{{-17,16,17,6,-9,-7,1,0,0,0},{-5,1,2,3,3,2,0,1,0,0},{1,5,-5,-1,0,-3,0,0,1,0},{4,-5,-5,-1,4,3,0,0,0,1}}
In[]:=
redmat=RowReduce[augmat];soln=FirstCase[redmat,{Repeated[0,{Length[bigvec]}],1,sol__}:>{sol}]
Out[]=
{-1,0,3}
Part 3: Find common solution mod integer n
Part 3: Find common solution mod integer n
We change the right hand side vectors so that we only have solutions mod 10.
In[]:=
SeedRandom[1234];newvec1=vec1+10*RandomInteger[{-5,5},3]newvec2=vec2+10*RandomInteger[{-5,5},3]mat1.solc-newvec1mat2.solc-newvec2
Out[]=
{-33,-6,23}
Out[]=
{4,59,-43}
Out[]=
{50,-10,-40}
Out[]=
{-10,-50,50}
Set up is similar to part 2 but now we augment bigmat below by 10 times an identity matrix of same number of columns, before augmenting by an identity matrix on the right. And because we must work over integers we use the Hermite Normal Form (HNF) instead of the Reduced Echelon Form. So we require HermiteDecomposition in place of RowReduce.
In[]:=
newbigvec=Join[newvec1,newvec2];newbigmat0=Join[-{newbigvec},bigmat0];newbigmat=Join[newbigmat0,10*IdentityMatrix[Length[newbigvec]]];newaugmat=Join[newbigmat,IdentityMatrix[Length[newbigmat]],2]
Out[]=
{{33,6,-23,-4,-59,43,1,0,0,0,0,0,0,0,0,0},{-5,1,2,3,3,2,0,1,0,0,0,0,0,0,0,0},{1,5,-5,-1,0,-3,0,0,1,0,0,0,0,0,0,0},{4,-5,-5,-1,4,3,0,0,0,1,0,0,0,0,0,0},{10,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0},{0,10,0,0,0,0,0,0,0,0,0,1,0,0,0,0},{0,0,10,0,0,0,0,0,0,0,0,0,1,0,0,0},{0,0,0,10,0,0,0,0,0,0,0,0,0,1,0,0},{0,0,0,0,10,0,0,0,0,0,0,0,0,0,1,0},{0,0,0,0,0,10,0,0,0,0,0,0,0,0,0,1}}
In[]:=
{uu,hnf}=HermiteDecomposition[newaugmat];solnmod=FirstCase[hnf,{Repeated[0,{Length[newbigvec]}],1,sol__}:>{sol}[[1;;Length[mat1[[1]]]]]]
Out[]=
{9,0,3}
Clearly this is a correct solution mod 10. We check though.
In[]:=
Mod[mat1.solnmod-newvec1,10]Mod[mat2.solnmod-newvec2,10]
Out[]=
{0,0,0}
Out[]=
{0,0,0}
Part 4: Extend to work over rationals
Part 4: Extend to work over rationals
This is the tricky part. We can clear denominators to work over the integers. For the example at hand we need to do this with the right hand sides. We use the least common multiple of the denominators (lcm) of course. We then augment below by that lcm times an appropriate identity matrix. If we get a row in the HNF matrix that has zeros in the slots where the right hand side vectors were, and a one in the next position, then subsequent elements in that row, divided by the lcm, will give the solution vector.
The catch is that we might not obtain such a vector in the HNF.
In[]:=
mats={{{-1,-1,0},{-31,-32,-32},{30,31,30}},{{-121,-109,-120},{-101,-91,-100},{210,189,208}}};vecs={{17/8,143/4,-44},{65/8,1289/8,137/8}};
In[]:=
bigmat0=Apply[Join[##,2]&,Map[Transpose,mats]];bigvec=-Flatten[vecs];lcm=LCM@@Denominator[bigvec]bigmat1=Join[lcm*{bigvec},bigmat0];bigmat2=Join[bigmat1,lcm*IdentityMatrix[Length[bigmat1[[1]]]]];bigmat=Join[bigmat2,IdentityMatrix[Length[bigmat2]],2];{uu,hh}=HermiteDecomposition[bigmat];keeprow=FirstCase[hh,{Repeated[0,{Length[bigvec]}],nz:Except[0],rest__}:>{nz,rest}]
Out[]=
8
Out[]=
{2,4,2,3,5,107,-122,149,433,-196}
What this shows is that, in order to annihilate the lcm times the combined right hand side vectors modulo the lcm (8), we have to multiply those vectors by 2. Which means we will not have a solution modulo the integers, but only modulo the half-integers.
What this means is that we will require a denominator that is a nontrivial multiple of the lcm, and moreover we know that multiplier is 2. So we repeat, now multiplying right hand side vectors and lower augmentation matrix by 2 times the lcm rather than just the lcm. Programmatically it is done as below.
In[]:=
mult=First[keeprow]newmult=mult*lcm;bigmat1new=Join[newmult*{bigvec},bigmat0];bigmat2new=Join[bigmat1new,newmult*IdentityMatrix[Length[bigmat1new[[1]]]]];bigmatnew=Join[bigmat2new,IdentityMatrix[Length[bigmat2new]],2];{uunew,hhnew}=HermiteDecomposition[bigmatnew];solnrow=FirstCase[hhnew,{Repeated[0,{Length[bigvec]}],nz:Except[0],rest__}:>{nz,rest}[[1;;Length[mats[[1,1]]]+1]]]
Out[]=
2
Out[]=
{1,12,2,3}
This means that we require {12,2,3} times the rows of the transposed, joined matrices to annihilate 2 times the lcm times the joined right-hand-side vector modulo 2 times the lcm. Which means our solution, mod 1, is {12,2,3}/16.
In[]:=
mod1soln=solnrow[[2;;]]/(2*lcm)mats[[1]].mod1soln-vecs[[1]]mats[[2]].mod1soln-vecs[[2]]
Out[]=
,,
3
4
1
8
3
16
Out[]=
{-3,-69,76}
Out[]=
{-135,-267,203}
Since these last are both integer vectors we have found zeros mod 1, so the claimed common solution is verified as correct.