3

我正在尝试求解方程

x1 + x2 + x3 + .... + xn = 1

其中 all 的值xi被限制为[0, 0.1, 0.2, ..., 0.9, 1]。

目前，我正在通过首先生成一个 n 维数组来解决这个问题mat，其中每个元素位置的值是轴值的总和，其变化范围为axisValues = 0:0.1:1：

mat(i,j,k,...,q) = axisValues(i) + axisValues(j) + ... + axisValues(q).

然后我搜索结果数组中等于一的所有条目。代码（如下所示以进一步说明）运行良好，并且已针对多达 5 个维度进行了测试。问题是，运行时间呈指数增长，我需要脚本在多个维度上工作。

clear all
dim = 2; % The dimension of the matrix is defined here. The script has been tested for dim ≤ 5
fractiles(:,1) = [0:0.1:1]; % Produces a vector containing the initial axis elements, which will be used to calculate the matrix elements
fractiles = repmat(fractiles,1,dim); % multiplies the vector to supply dim rows with the axis elements 0:0.1:1. These elements will be changed later, but the symmetry will remain the same.
dim_len = repmat(size(fractiles,1),1,size(fractiles,2)); % Here, the length of the dimensions is checked, which his needed to initialize the matrix mat, which will be filled with the axis element sums
mat = zeros(dim_len); % Here the matrix mat is initialized
Sub=cell(1,dim);
mat_size = size(mat);
% The following for loop fills each matrix elements of the dim dimensional matrix mat with the sum of the corresponding dim axis elements.
for ind=1:numel(mat)
    [Sub{:}]=ind2sub(mat_size,ind);
    SubMatrix=cell2mat(Sub);
    sum_indices = 0;
    for j = 1:dim
        sum_indices = sum_indices+fractiles(SubMatrix(j),j);
    end
    mat(ind) = sum_indices;
end
Ind_ones = find(mat==1); % Finally, the matrix elements equal to one are found.

我觉得以下使用问题对称性的想法可能有助于显着减少计算时间：

对于二维矩阵，所有满足上述条件的条目都位于从mat(1,11)到的对角线上mat(11,1)，即从 的最大值x1到 的最大值x2。

对于 3D 矩阵，所有条目都满足位于通过mat(1,1,11)、mat(1,11,1)、的对角线平面上的条件mat(11,1,1)，即从 的最大值x1到x2的最大值x3。

对于更高维也是如此：所有感兴趣的矩阵元素都位于一个n-1维超平面上，该维超平面固定在每个维中的最高轴值上。

问题是：有没有办法直接确定这些n-1维度超平面上元素的索引？如果是这样，整个问题可以一步解决，而无需计算 n 维矩阵的所有条目，然后搜索感兴趣的条目。

4

1 回答 1

3

数学：

我们不采用超立方体方式，而是求解方程

x(1) + x(2) + ... + x(n) = 1

每个都x(i)可以改变[0, 1/k, 2/k, ... (k-1)/k, 1]。在您的情况下k将是 10，因为这将导致 percents [0, 10, 20, ... 90, 100]。乘以k这个对应于丢番图方程

x(1) + x(2) + ... + x(n) = k,

所有的x(i)变化在哪里[0, 1, 2, ... k-1, k]。

我们可以在 this 和带有重复的组合的组合概念之间建立一个双射。

维基百科文章甚至隐含地提到了声明的潜在双射：

大小为 k 的多重子集的数量就是丢番图方程 的非负整数解的数量x1 + x2 + ... + xn = k。

对于一个较小的示例，假设我们将使用k=3百分比来[0, 33, 66, 100]代替。给定集合的所有 k-多重组合{1,2,3}：

RepCombs = 
     1     1     1
     1     1     2
     1     1     3
     1     2     2
     1     2     3
     1     3     3
     2     2     2
     2     2     3
     2     3     3
     3     3     3

然后我们使用以下规则将这些映射到您的百分比：对于每一行i，如果条目是j，则将1/3百分比添加到相应的矩阵条目中M(i,j)。第一行将对应于[1/3 + 1/3 + 1/3, 0, 0] = [1,0,0]。此过程生成的整体矩阵如下所示：

M =
    1.0000         0         0
    0.6667    0.3333         0
    0.6667         0    0.3333
    0.3333    0.6667         0
    0.3333    0.3333    0.3333
    0.3333         0    0.6667
         0    1.0000         0
         0    0.6667    0.3333
         0    0.3333    0.6667
         0         0    1.0000

代码：

现在对于生成所有这些的 MATLAB 代码：我使用这个答案中的函数nmultichoosek并实现我们accumarray的目标：

function M = possibleMixturesOfNSubstances(N, percentageSteps)
RepCombs = nmultichoosek(1:N, percentageSteps);
numCombs = size(RepCombs,1);
M = accumarray([repmat((1:numCombs).', percentageSteps, 1), RepCombs(:)], 1/percentageSteps, [numCombs, N]);

如果您想要百分比[0, 10, ... 90, 100]并且有 4 种物质，请使用以下方法调用此函数possibleMixturesOfNSubstances(4,10)

于 2015-02-14T11:42:50.020 回答