MATLAB Central Newsreader - Combinatorics question

Combinatorics question

Greg von Winckel — Thu, 29 Nov 2007 16:37:28 -0500

What would be an efficient way to compute all vectors with
natural number elements of length n<N such that the l^1 norm
of the vector is N?

Re: Combinatorics question

roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) — Thu, 29 Nov 2007 17:24:49 -0500

In article <fimps8$7in$1@fred.mathworks.com>,
Greg von Winckel <gregvw@gmail.com> wrote:
>What would be an efficient way to compute all vectors with
>natural number elements of length n<N such that the l^1 norm
>of the vector is N?

You mean "the knapsack problem" with a fixed number of weights
(but the size of each weight is limited only to natural numbers
1 thru N-n+1) ?
--
"All is vanity." -- Ecclesiastes

Re: Combinatorics question

John D'Errico — Thu, 29 Nov 2007 17:53:55 -0500

"Greg von Winckel" <gregvw@gmail.com> wrote in message
<fimps8$7in$1@fred.mathworks.com>...
> What would be an efficient way to compute all vectors with
> natural number elements of length n<N such that the l^1 norm
> of the vector is N?

My partitions code from the FEX is close to
what you are asking for. Is it efficient?
Reasonably so.

http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?
objectId=12009&objectType=FILE

The only flaw is that its set up to limit only the
quantity of any single number, not to fix the
total number of elements.

However, it should be reasonably efficient to
write a recursive code to solve your problem
using much the same scheme found in
partitions. For example:

% =============================
function list = parts(Nsum,Nel,MaximumElement)
% parts: find all decreasing positive integer vectors of length Nel that sum to
Nsum
% list = parts(Nsum,Nel)
%
% arguments:
% Nsum - scalar integer - defines the sum
%
% Nel - scalar integer - number of elements in the vector

if (nargin<3) || isempty(MaximumElement)
  MaximumElement = Nsum;
end

if (Nel == Nsum)
  % the simple case, where all elements are equal and unity
  list = ones(1,Nel);
  return
elseif (Nsum <= 0) || (Nel == 0)
  % either of these events cannot be satisfied by any list
  % of positive integers
  list = [];
  return
elseif (Nel == 1)
  % Only one element in the list - it must be Nsum
  if (MaximumElement >= Nsum)
    list = Nsum;
    return
  else
    % but that possibility is disallowed
    list = [];
    return
  end
elseif (Nel>Nsum)
  % another non-viable possibility
  list = [];
  return
end

% if we drop down to here, there are at least two elements, and
% at least one solution to be found.
list = cell(0,1);
for maxel = min(MaximumElement,Nsum - Nel + 1):-1:1
  % Assume that the first element in the list is maxel
  sublist = parts(Nsum - maxel,Nel-1,maxel);
  if ~isempty(sublist)
    r = size(sublist,1);
    list{end+1,1} = [repmat(maxel,r,1),sublist];
  end
end

% aggregate all solutions into a flat array
list = cell2mat(list);
% =============================

parts(10,4)
ans =
     7 1 1 1
     6 2 1 1
     5 3 1 1
     5 2 2 1
     4 4 1 1
     4 3 2 1
     4 2 2 2
     3 3 3 1
     3 3 2 2

parts(20,4)
ans =
    17 1 1 1
    16 2 1 1
    15 3 1 1
    15 2 2 1
    14 4 1 1
    14 3 2 1
    14 2 2 2
    13 5 1 1
    13 4 2 1
    13 3 3 1
    13 3 2 2
    12 6 1 1
    12 5 2 1
    12 4 3 1
    12 4 2 2
    12 3 3 2
    11 7 1 1
    11 6 2 1
    11 5 3 1
    11 5 2 2
    11 4 4 1
    11 4 3 2
    11 3 3 3
    10 8 1 1
    10 7 2 1
    10 6 3 1
    10 6 2 2
    10 5 4 1
    10 5 3 2
    10 4 4 2
    10 4 3 3
     9 9 1 1
     9 8 2 1
     9 7 3 1
     9 7 2 2
     9 6 4 1
     9 6 3 2
     9 5 5 1
     9 5 4 2
     9 5 3 3
     9 4 4 3
     8 8 3 1
     8 8 2 2
     8 7 4 1
     8 7 3 2
     8 6 5 1
     8 6 4 2
     8 6 3 3
     8 5 5 2
     8 5 4 3
     8 4 4 4
     7 7 5 1
     7 7 4 2
     7 7 3 3
     7 6 6 1
     7 6 5 2
     7 6 4 3
     7 5 5 3
     7 5 4 4
     6 6 6 2
     6 6 5 3
     6 6 4 4
     6 5 5 4
     5 5 5 5

Now I'd write a permutation code for each
row of this result. Expect this all to get nasty
exponentially fast for larger sums. parts(50,10)
found 16928 solutions, each of which can occur
in a multitude of sequences.

John

Re: Combinatorics question

Bruno Luong — Thu, 29 Nov 2007 18:37:03 -0500

"John D'Errico" <woodchips@rochester.rr.com> wrote in
> parts(10,4)
> ans =
> 7 1 1 1
> 6 2 1 1
> 5 3 1 1
> 5 2 2 1
> 4 4 1 1
> 4 3 2 1
> 4 2 2 2
> 3 3 3 1
> 3 3 2 2

Is it all combination? I find 84 possible solutions:

1 1 1 7
1 1 2 6
1 1 3 5
1 1 4 4
1 1 5 3
1 1 6 2
1 1 7 1
1 2 1 6
1 2 2 5
1 2 3 4
1 2 4 3
1 2 5 2
1 2 6 1
1 3 1 5
1 3 2 4
1 3 3 3
1 3 4 2
1 3 5 1
1 4 1 4
1 4 2 3
1 4 3 2
1 4 4 1
1 5 1 3
1 5 2 2
1 5 3 1
1 6 1 2
1 6 2 1
1 7 1 1
2 1 1 6
2 1 2 5
2 1 3 4
2 1 4 3
2 1 5 2
2 1 6 1
2 2 1 5
2 2 2 4
2 2 3 3
2 2 4 2
2 2 5 1
2 3 1 4
2 3 2 3
2 3 3 2
2 3 4 1
2 4 1 3
2 4 2 2
2 4 3 1
2 5 1 2
2 5 2 1
2 6 1 1
3 1 1 5
3 1 2 4
3 1 3 3
3 1 4 2
3 1 5 1
3 2 1 4
3 2 2 3
3 2 3 2
3 2 4 1
3 3 1 3
3 3 2 2
3 3 3 1
3 4 1 2
3 4 2 1
3 5 1 1
4 1 1 4
4 1 2 3
4 1 3 2
4 1 4 1
4 2 1 3
4 2 2 2
4 2 3 1
4 3 1 2
4 3 2 1
4 4 1 1
5 1 1 3
5 1 2 2
5 1 3 1
5 2 1 2
5 2 2 1
5 3 1 1
6 1 1 2
6 1 2 1
6 2 1 1
7 1 1 1

Bruno

Re: Combinatorics question

Bruno Luong — Thu, 29 Nov 2007 18:42:15 -0500

If you consider "natural number" starts from 0, call:
v=allVL1(n, L1);

If you consider "natural number" starts from 1, then call:
v=allVL1(n, L1-n)+1;

Bruno

function v=allVL1(n, L1)
% function v=allVL1(n, L1)
% INPUT
% n: length of the vector
% L1: desired L1 norm
% OUTPUT:
% v: m x n array such as sum(v,2)=L1
% all elements of v is naturel numbers {0,1,...}
% v contains all possible combination

if n==1
    v = L1;
else
    v1 = (0:L1)';
    vrest = arrayfun(@(j) allVL1(n-1, L1-j), v1, ...
        'UniformOutput', false);
    v = arrayfun(@(i) ...
        [repmat(v1(i),size(vrest{i},1),1) vrest{i}], ...
        (1:length(v1))',...
        'UniformOutput', false);
    v = cell2mat(v);
end

Re: Combinatorics question

Steven Lord — Thu, 29 Nov 2007 20:53:01 -0500

"Bruno Luong" <b.luong@fogale.fr> wrote in message
news:fin0sf$g6s$1@fred.mathworks.com...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in
>> parts(10,4)
>> ans =
>> 7 1 1 1
>> 6 2 1 1
>> 5 3 1 1
>> 5 2 2 1
>> 4 4 1 1
>> 4 3 2 1
>> 4 2 2 2
>> 3 3 3 1
>> 3 3 2 2
>
> Is it all combination? I find 84 possible solutions:
>
> 1 1 1 7
> 1 1 2 6
> 1 1 3 5
> 1 1 4 4
> 1 1 5 3
> 1 1 6 2
> 1 1 7 1
> 1 2 1 6

*snip*

Note that some of the elements of John's list appear multiple times in your
list, just with the elements permuted. For instance, [1 1 1 7] and [1 1 7
1] are each permutations of [7 1 1 1] in John's list. The OP's question
didn't specify whether the order of the elements in the vector was important
for his application.

--
Steve Lord
slord@mathworks.com

Re: Combinatorics question

Bruno Luong — Thu, 29 Nov 2007 21:41:05 -0500

Same function, slightly improved.

Bruno

function [v m]=allVL1(n, L1)
% function v=allVL1(n, L1)
% INPUT
% n: length of the vector
% L1: desired L1 norm
% OUTPUT:
% v: (m x n) array such as sum(v,2)=L1
% all elements of v is naturel numbers {0,1,...}
% v contains all possible combinations
% Remark:
% allVL1(n,L1-n)+1 for natural numbers defined as {1,2,...}

if n==1
    v = L1;
    m = 1;
else
    v1 = (0:L1)';
    [vrest m]= arrayfun(@(j) allVL1(n-1, L1-j), v1, ...
                        'UniformOutput', false);
    v1 = arrayfun(@(n,m) n+zeros(m,1), v1, [m{:}]', ...
                  'UniformOutput', false);
    v = [cell2mat(v1) cell2mat(vrest)];
    m = size(v,1);
end

Re: Combinatorics question

Bruno Luong — Thu, 29 Nov 2007 22:50:12 -0500

Probably the best I can do: ;-)

Bruno

function v=allVL1(n, L1, head)
% function v=allVL1(n, L1);
% INPUT
% n: length of the vector
% L1: desired L1 norm
% OUTPUT:
% v: (m x n) array such as sum(v,2)=L1
% all elements of v is naturel numbers {0,1,...}
% v contains all (=m) possible combinations
% Remark:
% allVL1(n,L1-n)+1 for natural numbers defined as {1,2,...}

if n==1
    v = L1;
else
    v = cell2mat(arrayfun(@(j) allVL1(n-1, L1-j, j), ...
                          (0:L1)', 'UniformOutput', false));
end

if nargin>=3
    v = [head+zeros(size(v,1),1) v];
end

Re: Combinatorics question

Greg von Winckel — Fri, 30 Nov 2007 10:10:50 -0500

Thanks for all the ideas. I will check out your codes. I do
indeed need every permutation of the solution vectors.

I'm working on some sparse grid codes in MATLAB and the
interpolants on a sparse grid of order N are constructed from
all possible combinations of tensor product interpolants
where the indices of the constituent dimensions add up to N.

Thanks for the suggestions!

Greg

Re: Combinatorics question

Bruno Luong — Fri, 30 Nov 2007 10:43:54 -0500

I have submitted my code in FEX. You can also select the
criteria L1 norm by ==, or <= or < to some quantity.

Bruno

Re: Combinatorics question

Greg von Winckel — Fri, 30 Nov 2007 10:57:23 -0500

Bruno, your code is very compact and looks good, although I
should mention that zero is not a natural number. I can make
the necessary modifications though.

"Bruno Luong" <b.luong@fogale.fr> wrote in message
<fiopha$jb7$1@fred.mathworks.com>...
> I have submitted my code in FEX. You can also select the
> criteria L1 norm by ==, or <= or < to some quantity.
>
> Bruno

Re: Combinatorics question

Bruno Luong — Fri, 30 Nov 2007 11:12:35 -0500

"Greg von Winckel" <gregvw@gmail.com> wrote in message
<fioqaj$rrr$1@fred.mathworks.com>...
> Bruno, your code is very compact and looks good, although I
> should mention that zero is not a natural number. I can make
> the necessary modifications though.
>

Sure, feel free to modify it. Or as I write above, you could
also call the original function:

V = allVL1(n,L1-n)+1

and this will return the vectors of natural numbers starting
from 1, and satified sum(V,2)=L1.

Bruno