Thread Subject: Integer solutions for a_1*x_1+a_2*x_2+a_3*x_3+ ...+a_4*x_n = k

Hello,

I am wondering how to write a function for solving

a_1*x_1+a_2*x_2+a_3*x_3+ ...+a_4*x_n = k

where
Inputs:
n, k, and vector of non-negaive integer coefficients a=[a_1, a_2, ...a_n].
k is some small positive integer, for example, k= 3, or k= 6, n is rather large positve integer, for example,, n = 20, or n= 60 or in some rare cases even n=100.

Output:
The unknowns x_1, x_2, ...x_n are non-negative integers.

I have asked a simplier version of this problem here, where a=[1, 1, ...1], and got lots of help from Walter, Bruno, Roger, John (thank you).

But this problem... seems really difficult, is it even possible to solve it in Matlab? Any ideas?

Thank you for your help!
Oriole






Thank you!!

"Oriole " <oriole_ni@hotmail.com> wrote in message <gfun9a$41t$1@fred.mathworks.com>...
> Hello,
>
> I am wondering how to write a function for solving
>
> a_1*x_1+a_2*x_2+a_3*x_3+ ...+a_4*x_n = k
>
> where
> Inputs:
> n, k, and vector of non-negaive integer coefficients a=[a_1, a_2, ...a_n].
> k is some small positive integer, for example, k= 3, or k= 6, n is rather large positve integer, for example,, n = 20, or n= 60 or in some rare cases even n=100.
>
> Output:
> The unknowns x_1, x_2, ...x_n are non-negative integers.
>
> I have asked a simplier version of this problem here, where a=[1, 1, ...1], and got lots of help from Walter, Bruno, Roger, John (thank you).
>
> But this problem... seems really difficult, is it even possible to solve it in Matlab? Any ideas?
>
> Thank you for your help!
> Oriole


As I'm sure you realize, this problem may not have a solution.

However, if it does have one, it looks like it could be solved efficiently by modeling the problem as a Markov Decision Process and using dynamic programming. Not sure if you have a background in that. If not, a useful reference is Chapt 4 in

"Markov Decision Processes: Discrete Stochastic Dynamic Programming" by M.L. Puterman

Basically, the idea is to model this as a sequence of n decisions in a controlled Markov chain as follows:

The initial state is s_0=k.

The 1st decision is the chosen value of x_1 and causes a state transition to s_1=k-a_1*x_1 with cost=0.

Thereafter, the m-th state transition is from s_(m-1) to s_m=s_(m-1)-a_m*k_m also with cost 0.

After n decisions you will have a terminal cost equal to s_n, which is the residual amount

a_1*x_1+a_2*x_2+a_3*x_3+ ...+a_4*x_n - k

that you want to minimize.

This terminal cost will also be your total cost, because we've set all other decision costs to zero. Dynamic programming will minimize this total cost and give you a solution x_1...x_n (if it exists).

Hope this helps.

"Oriole " <oriole_ni@hotmail.com> wrote in message <gfun9a$41t$1@fred.mathworks.com>...
> Hello,
>
> I am wondering how to write a function for solving
>
> a_1*x_1+a_2*x_2+a_3*x_3+ ...+a_4*x_n = k
>
> where
> Inputs:
> n, k, and vector of non-negaive integer coefficients a=[a_1, a_2, ...a_n].
> k is some small positive integer, for example, k= 3, or k= 6, n is rather large positve integer, for example,, n = 20, or n= 60 or in some rare cases even n=100.
>
> Output:
> The unknowns x_1, x_2, ...x_n are non-negative integers.
>
> I have asked a simplier version of this problem here, where a=[1, 1, ...1], and got lots of help from Walter, Bruno, Roger, John (thank you).
>
> But this problem... seems really difficult, is it even possible to solve it in Matlab? Any ideas?
>
> Thank you for your help!
> Oriole

This is a Diophantine equation, as is the partitions
problem, which is just a special case of the
general Diophantine equation.

For n very large, this is a nasty problem. You can still
solve it recursively as we have shown you, but that
may be computationally expensive.

John

> Thereafter, the m-th state transition is from s_(m-1) to s_m=s_(m-1)-a_m*k_m also with cost 0.

Sorry, this should be s_(m-1) to s_m=s_(m-1)-a_m*x_m

Note also that the decision variables will be restricted to 1<=x_m<=k

> After n decisions you will have a terminal cost equal to s_n, which is the residual amount
>
> a_1*x_1+a_2*x_2+a_3*x_3+ ...+a_4*x_n - k

And this should be

k-a_1*x_1+a_2*x_2+a_3*x_3+ ...+a_4*x_n

....But hopefully you get the idea

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <gfuq8v$n6o$1@fred.mathworks.com>...
> "Oriole " <oriole_ni@hotmail.com> wrote in message <gfun9a$41t$1@fred.mathworks.com>...
> > Hello,
> >
> > I am wondering how to write a function for solving
> >
> > a_1*x_1+a_2*x_2+a_3*x_3+ ...+a_4*x_n = k
> >
> > where
> > Inputs:
> > n, k, and vector of non-negaive integer coefficients a=[a_1, a_2, ...a_n].
> > k is some small positive integer, for example, k= 3, or k= 6, n is rather large positve integer, for example,, n = 20, or n= 60 or in some rare cases even n=100.
> >
> > Output:
> > The unknowns x_1, x_2, ...x_n are non-negative integers.
> >
> > I have asked a simplier version of this problem here, where a=[1, 1, ...1], and got lots of help from Walter, Bruno, Roger, John (thank you).
> >
> > But this problem... seems really difficult, is it even possible to solve it in Matlab? Any ideas?
> >
> > Thank you for your help!
> > Oriole
>
> This is a Diophantine equation, as is the partitions
> problem, which is just a special case of the
> general Diophantine equation.
>
> For n very large, this is a nasty problem. You can still
> solve it recursively as we have shown you, but that
> may be computationally expensive.
>
> John

For kicks, I wrote a simple recursive code to compute
the integer solutions. For a problem with 15 coefficients,
it runs reasonably quickly, for small k (k = 5 here.)

A = rem(1:6,3)+1
A =
     2 3 1 2 3 1


Each row of x is a solution.

>> x = diophantine(A,5,0:5)
x =
     1 1 0 0 0 0
     2 0 1 0 0 0
     0 1 2 0 0 0
     1 0 3 0 0 0
     0 0 5 0 0 0
     0 1 0 1 0 0
     1 0 1 1 0 0
     0 0 3 1 0 0
     0 0 1 2 0 0
     1 0 0 0 1 0
     0 0 2 0 1 0
     0 0 0 1 1 0
     2 0 0 0 0 1
     0 1 1 0 0 1
     1 0 2 0 0 1
     0 0 4 0 0 1
     1 0 0 1 0 1
     0 0 2 1 0 1
     0 0 0 2 0 1
     0 0 1 0 1 1
     0 1 0 0 0 2
     1 0 1 0 0 2
     0 0 3 0 0 2
     0 0 1 1 0 2
     0 0 0 0 1 2
     1 0 0 0 0 3
     0 0 2 0 0 3
     0 0 0 1 0 3
     0 0 1 0 0 4
     0 0 0 0 0 5

>> tic,x = diophantine(rem(1:15,3)+1,5,0:5);toc
Elapsed time is 1.112639 seconds.
>> size(x)
ans =
   476 15

tic,x = diophantine(rem(1:20,3)+1,5,0:5);toc
Elapsed time is 2.848407 seconds.
>> size(x)
ans =
        1008 20

>> tic,x = diophantine(rem(1:25,3)+1,5,0:5);toc
Elapsed time is 8.448537 seconds.
>> size(x)
ans =
        2592 25

>> tic,x = diophantine(rem(1:30,3)+1,5,0:5);toc
Elapsed time is 19.793219 seconds.
>> size(x)
ans =
        5402 30


Expect this to be slow for most serious problems.


function xsol = diophantine(A,c,xset)
% solves the diophantine equation dot(A,x) = c, where x is a member of xset
%
% arguments
% A - row vector of length p, integer coefficients
% c - scalar (integer) constant
% xset - set to sample from, sampling with replacement
%
% xsol - mxp array of solutions found

% check sizes, shapes
A = A(:)';
p = length(A);
xset = xset(:);
nx = length(xset);

xsol = zeros(0,p);
if length(A) == 1
  ind = find((A*xset - c) == 0);
  if ~isempty(ind)
    xsol = xset(ind);
  end
else
  % recursive calls
  for i = 1:nx
    if (c-A(end)*xset(i))>=0
      xsub = diophantine(A(1:(end-1)),c-A(end)*xset(i),xset);
      nsub = size(xsub,1);
      xsol = [xsol;[xsub,repmat(xset(i),nsub,1)]];
    end
  end
end

Have fun,
John

"Oriole " <oriole_ni@hotmail.com> wrote in message <gfun9a$41t$1@fred.mathworks.com>...
> I am wondering how to write a function for solving
> ........
> Oriole

  It looks as though John has already done the problem for you. It is the way I would have done it, using recursion. For the small k you mention it should not lead to excessive computation times.

  About the only thing original I can contribute is that your restriction of the a values to being only non-negative, as opposed to positive, integers seems inappropriate. Any time you have a solution with one of the a's, say a_i, equal to zero, there will be infinitely many solutions with all possible non-negative values for x_i and that would lead to a failure of John's recursive algorithm. You should really restrict the a's to being just the positive integers.

Roger Stafford

"Oriole " <oriole_ni@hotmail.com> wrote in message <gfun9a$41t$1@fred.mathworks.com>...
> Hello,
>
> I am wondering how to write a function for solving
>
> a_1*x_1+a_2*x_2+a_3*x_3+ ...+a_4*x_n = k
>
> where
> Inputs:
> n, k, and vector of non-negaive integer coefficients a=[a_1, a_2, ...a_n].
> k is some small positive integer, for example, k= 3, or k= 6, n is rather large positve integer, for example,, n = 20, or n= 60 or in some rare cases even n=100.
>
> Output:
> The unknowns x_1, x_2, ...x_n are non-negative integers.

It wasn't clear to me whether you were looking for one solution or all of them.
If you require only one solution, my dynamic programming idea is applicable and I've implemented it below. It is quite fast, even for n=100


>> aa=round(rand(1,100)*10+2);

>> tic; xx=dioph(aa,50); toc %test for n=100 and k=50
Elapsed time is 0.000800 seconds.


Of course, I don't know if my aa test vector is representative of your application...








###############################

function xx=dioph(aa,kk)
%Finds non-negative integer solution xx to the equation
%
% aa'*xx=kk
%
% where aa is a strictly positive integer vector and kk is a non-negative
% scalar

if any(aa-round(aa))
    error 'Not integer valued a'
end

if any(aa<=0)
    error 'Should be strictly positive a'
end

nn=length(aa);
states=(-1:kk).';
Nstates=kk+2;
Nactions=kk+1;

Smap=repmat(states,[1,Nactions]);
Dmap=repmat((0:kk),[Nstates,1]);

TerminalCost=states; TerminalCost(1)=2*kk;

Value=TerminalCost;
StateAction=zeros(Nstates,Nactions);

for ii=1:nn
        
    
    M=Smap-aa(ii)*Dmap;
    M(M<0)=-1;
    M=M+2;
    
 
    [Value,StateAction(:,ii)]=min(Value(M),[],2);
    StateAction(:,ii)=StateAction(:,ii)-1;
    
    
   
    if (Value(end)==0)
        LastStage=ii;
       break;
    end
    
end

xx=zeros(1,nn);
TheState=kk;

for jj=LastStage:-1:1

 TheAction=StateAction(TheState+2,jj);
 xx(jj)=TheAction;
 TheState=TheState-aa(jj)*TheAction;
 
    
end

Hi Oriole,

I have wrote a code for this problem few months backs. Here is the main function. Please save it in a file intlin.m. This is recursion algo, but call for integer common division properties and linear programming to accelerate the computation.

A test file will follow.

Bruno

PS: this is from an old thread
http://www.mathworks.com/matlabcentral/newsreader/view_thread/172882

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function x = intlin(a, b)
% function x = intlin(a, b)
%
% PURPOSE: solving integer equation:
%
% x(1)*a(1) + x(2)*a(2) + ... + x(n)*a(n) = b
% x >=0
%
% Input: a (1 x n) integer (could be negative)
% b (1 x 1) integer (could be negative)
% Output: x (n x 1) integer such that
% a*x = b, x>=0
%
% Author: Bruno Luong
% Last update: 24/July/2008

% Default value, empty solution
x = zeros(0,length(a));

if ~isscalar(b)
    fprintf('b must be scalar\n');
    return
end

% cast to double, and reshape as long thin array
a = double(a(:));
b = double(b);

if ~all(mod(a,1)==0) || mod(b,1)~=0
    fprintf('a and b must be integers\n');
    return
end

% d*a' = g
[g d] = gcdn(a);

if mod(b,g)==0
    an = a/g;
    bn = b/g;
    
    % (d*a' = b) OR (d*an' = bn)
    d = bn*d;
    
    %
    % General gcd final solution would be
    % x = d + k;
    % with k such that: k*an' = 0, and
    % k>=-d (<=> x>=0)
    kmin = ceil(-d);
    kmax = +inf(size(a));

    % Find all k such that k.an = 0
    k = allintlin0(an, kmin, kmax);
    
    x = bsxfun(@plus, d, k);
    
else % mode(b,g)~=0
    fprintf('WARNING: there is no solution\n');
end

end


function x = allintlin0(a, lower, upper)
%
% x, a, lower, upper are n-dimensional vector
% List all interger x such that
% x.a = 0
% lower <= x <= upper
%

a = a(:);
n = length(a);

lower = lower(:);
upper = upper(:);
% Adjust upper and lower bounds by LP
L = lower;
U = upper;
b = 0;
epsilon = 1e-6;
for k=1:n
    cost = basis(k,n);
    % Beware, this is Bruno's linprog, not MATLAB one in optimization
    % tool box, the result may be different.
    sol = linprog(cost, [], [], a', b, L, U);
    if ~all(isinf(sol))
        L(k) = max(L(k),sol(k)-epsilon);
    end
    cost = -basis(k,n);
    sol = linprog(cost, [], [], a', b, L, U);
    if ~all(isinf(sol))
        U(k) = min(U(k),sol(k)+epsilon);
    end
end
L = ceil(L);
U = floor(U);

if all(~isinf(L)) && all(~isinf(U))
    maxcount = Inf;
else % Limit the number of solutions that will be listed
    % NOTE: set maxcount to finite value doesn't work as efficienly,
    % as expected because the recursive engine might spend much of
    % CPU time to look for unvalid solutions

    % maxcount = 100;
    fprintf('There is infinity of solutions\n');
    x = NaN;
    return
end

% Call the engine
count = 0;
x = ilinengine(a, L, U, b, count, maxcount);
if maxcount<inf
    fprintf('WARNING: only %d solutions will be provided\n', maxcount);
    x = x(1:min(maxcount,end),:); % clipping
end

end

function v = basis(k,n) % generate a k-th basis vector of dimension n
v = zeros(n,1);
v(k) = 1;
end

% Solver engine for integer x
% a*x = b
% lower <= x <= upper
% RESTRICTION:
% a must be primary array, i.e., they greatest common divisor is one
function x = ilinengine(a, lower, upper, b, count, maxcount)

% default value, empty result
x = zeros(0,length(a));

if count>maxcount
    return
end

% Trivial solution with one variable
if length(a)==1
    if mod(b,a(1))==0
        xtmp = b / a(1);
        if (xtmp >= lower) && (xtmp <=upper)
            x = xtmp;
        end
    end
    return
end

% preliminary check where as the sum b is possible
% This check is to speed up (?), and deos not affect the result
as = bsxfun(@times,a,[lower upper]);
as = sum(sort(as,2),1);
if b<as(1) || b>as(2)
    return
end

% Greatest common divisor for the tail
g = gcdn(a(2:end));
% find r1, the inverse of a(1) in g-modulo group
[uno r1] = gcd(a(1),g);
clear uno; % should be 1

r=r1*mod(b,g);
s = (b - a(1)*r)/g;
a(2:end) = a(2:end)/g; % the tail is primary among them

kmin = ceil((lower(1)-r)/g);
kmax = floor((upper(1)-r)/g);

    % perform a basic step
    function newcount = kstep(k)
        % Recursive call
        xk = ilinengine(a(2:end), lower(2:end), upper(2:end), ...
                        s-k*a(1), count, maxcount);
        nxk = size(xk,1);
        x = [x; ...
            (r+g*k)+zeros(nxk,1) xk]; % append the new solutions
        newcount = count + nxk; % adjust the counter
    end

if isinf(kmin) && isinf(kmax) % k is unbounded
    for absk=0:inf
        for k=unique([-absk absk])
            if kstep(k)>maxcount
                break
            end
        end
        if count>maxcount
            break
        end
    end
elseif isinf(kmin) % k has upper bound, but no lower bound
    for k=kmax:-1:-inf
        if kstep(k)>maxcount
            break
        end
    end
else % k has lower bound
    for k=kmin:kmax
        if kstep(k)>maxcount
            break
        end
    end
end

end

function [g varargout] = gcdn(varargin)
% function g = gcdn(a1, a2, ..., an);
%
% Return g, Greatest common divisor of a1, ... an
%
% [g c1 c2 ... cn]=gcdn(a1, a2, ..., an)
% return also c1, ..., cn
% So that a1*c1 + ... an*cn = g
%
% Compact calling form:
% [g c]=gcdn(a1, a2, ..., an) or [g c]=gcdn(a)
% assumes a and c are array
%

if nargin<2
    a = varargin{1};
    if length(a)<2
        if a
            g = abs(a);
            if nargout>=2
                varargout{1}=sign(a);
            end
            return
        else
           error('gcdn cannot compute for a = 0');
        end
    end
    a = reshape(a,1,[]);
else
    % Put all numbers in array
    a=cell2mat(varargin);
end

g=a(1);
c=zeros(size(a));
c(1)=1;
for k=2:length(a)
    [g cg c(k)]= gcd(g, a(k));
    c(1:k-1) = c(1:k-1)*cg;
end

if nargout>=2
    switch (nargout-1)
        case 1,
            varargout{1}=c;
        case length(c)
            varargout=num2cell(c);
        otherwise
             error('number of output is incompatible with input');
    end
end

end

% Test program:

% This example would take more or less a minute to solve
a=[770 105 60 85 50];

b = 9435;

%
% solving
% c1*a1 + ... cn*an = b
% c >=0

c = intlin(a, b);

if isnan(c)
elseif ~isempty(c)

    fprintf('\nSolutions c:\n\n');
    disp(c);
    fprintf('\nNumber of solutions = %d\n\n', size(c,1));
    
    if exist('cc','var') && ~ismember(cc,c,'rows')
        fprintf('Solver miss the initial solution\n');
        %save test_inlin_debug.mat a cc b
    end
   
    bverif = c*a';
    bverifu = unique(bverif);
    
    if length(bverifu)~=1 || bverifu~=b
        fprintf('Wrong solution(s)\n');
        disp(c(bverif~=b,:));
        %save test_inlin_debug.mat a cc b
    end
    
else % mode(b,g)==0
    fprintf('There is no solution\n');
end

% Bruno