最小化Matlab中数组列的总和
Minimising the sum of array columns in Matlab
我有一个大数组（大约 250,000 x 10）。每行包含 1 或 -1。例如：
data(1, :) = [1, -1, -1, -1, -1, -1, -1, -1, 1, -1];
我需要 select 组 n 行，这样 列的绝对总和的平均值 被最小化（尽可能接近零）。因此，在这个玩具示例中，其中 n=2：
[ 1 1 1 1]
[-1 -1 -1 -1]
[-1 1 -1 1]
我会 select 第 1 行和第 2 行，因为它们总和为 [0 0 0 0]（平均值 0），这是 n=2 时可能的最小值。
我尝试了下面建议的方法（寻找互补对），但对于我的数据集，这只能形成 23k 行的平衡子集。因此，我需要一个近似值，它生成大小为 n 行的子集，但具有列的绝对和的最小均值。
到目前为止我发现的最佳方法如下：选择一个起始子集，将剩余的每一行迭代添加到基数，如果它提高了列的绝对和的平均值则保留它。这是非常粗糙的，我相信有更好的方法。它也很容易卡在错误的最小值，所以需要添加一个意外事件：
shuffle = randperm(size(data));
data_shuffled = data(shuffle, :);
base = data_shuffled(1:30000, :);
pool = data_shuffled(30001:end, :);
best_mean = mean(abs(sum(base, 1)));
best_matrix = base;
n = 100000;
for k = 1:20
for i = 1:size(pool, 1)
temp = pool (i, :);
if(~isnan(temp(1)))
temp_sum = sum(temp, 1);
new_sum = temp_sum + sum(best, 1);
temp_mean = mean(abs(new_sum));
if(temp_mean < best_mean)
best_mean = temp_mean;
best_matrix = vertcat(best_matrix, temp);
pool(i, :) = NaN(1, 10);
end
end
end
if(size(best_matrix, 1) > n)
return
end
end
这实现了 ~17,000 列的绝对总和的平均值，这还不错。重复使用不同的种子可能会有所改善。
理想情况下，我不会将新元素添加到 best_matrix 的末尾，而是 将它 替换为一些元素，以实现最佳改进。
更新：我不想给出数据集的具体细节，因为所有解决方案都应该适用于指定格式的任何矩阵。
感谢大家的贡献！
不幸的是，这个问题超出了常规（连续）优化的范围。你的问题，可以参数化如下：
min_{S∈S_n} Σ_{j∈S}|Σ_i data_ji|
其中S_n是索引j∈{0,...,250000}的n个元素组合的集合，也可以重写为非常类似x中的正则二次整数规划问题：
min_x x'* data *data' *x
0<=x<=1 and x*1=n
其中 data 是您的 250000*10 矩阵，x 是我们正在寻找的 250000*1 组合向量。 （现在我们优化的是平方和，而不是绝对值的和...）
这道题被证明是NP-hard，也就是说要求全局最小值，必须遍历n种所有可能的组合，抽取250000种可能，这等于二项式系数（250000 n），等于250000!/(n!*(250000-n)!)...
祝你好运...;)
编辑
如果您要启发式地解决这个问题，因为我想您将需要一个解决方案，请使用启发式 here 而不是您的方法。
正如其他人所说，最佳解决方案可能是不可能的，所以我将专注于特定情况。
首先我假设每列的分布是独立的。
然后我在累加器 space 上工作，以减少数据大小并加快代码速度。
我通过将每个 -1 视为 0 并将每一行视为一个数字，并添加 1 以避免使用 0 作为索引，例如：
data(1,:)=[-1 1 -1 1 -1 1 -1 1 -1 1] -> '0101010101' -> 341 -> 342
有了这个我们可以累积数据为：
function accum=mat2accum(data)
[~,n]=size(data);
indexes=bin2dec(num2str((data+1)/2))+1;
accum=accumarray(indexes,1,[2^n 1]);
我考虑的第一种情况是当每列的总和与数据的大小相比较小时，这意味着所有列中的 1 和 -1 的数量相似。
sum(data) << size(data)
对于这种情况，您可以找到所有相互抵消的对，例如：
data(1,:)=[-1 1 -1 1 -1 1 -1 1 -1 1] -> '0101010101' -> 341 -> 342
data(2,:)=[1 -1 1 -1 1 -1 1 -1 1 -1] -> '1010101010' -> 682 -> 683
并且我们知道每一对都将位于累加器索引中的镜像位置，因此我们可以通过以下方式获得所有可能的对：
function [accumpairs, accumleft]=getpairs(accum)
accumpairs=min([accum,accum(end:-1:1)],[],2);
accumleft=accum-accumpairs;
使用随机生成的数据，我能够在一组 250k 行中获得 >100k 对，并且每对的子集在每列中的总和等于零。因此，如果 1 和 -1 均匀分布，这可能就足够了。
我考虑的第二种情况是每列的总和远不为零，这意味着 1 和 -1 之间存在很大的不成比例。
abs(sum(data)) >> 0
通过反转总和为负的每一列，这不会影响数据，因为最后可以再次反转这些列。可以强制不成比例为 1 多于 -1。并且通过首先提取该数据的可能对，不成比例更加明显。
有了这样准备的数据，就可以将问题视为最小化所需集合中 1 的数量。为此，我们首先将可能的指标随机化，然后计算并排序每个指标的汉明权重（二进制表示中1的个数），然后收集汉明权重尽可能小的数据。
function [accumlast,accumleft]=resto(accum,m)
[N,~]=size(accum);
columns=log2(N);
indexes=randperm(N)'; %'
[~,I]=sort(sum((double(dec2bin(indexes-1,columns))-48),2));
accumlast=zeros(N,1);
for k=indexes(I)' %'
accumlast(k)=accum(k);
if sum(accumlast)>=m
break
end
end
accumleft=accum-accumlast;
使用随机生成的数据，其中 1 的数量大约是 -1 的 2 倍，并且每列的总和约为 80k，我可以找到 100k 数据的子集，每列的总和约为 5k。
第三种情况，有些列的总和接近于零，有些则不是。在这种情况下，您将列分为具有大和的列和具有小和的列，然后按大和列的汉明权重对数据进行排序，并在每个大列索引中获取小和列对.这将为大和列的每个索引创建一个矩阵，其中包含可能对的数量、不成对的行数以及小列的不成对行的总和。
现在您可以使用该信息来保存 运行 总和，并查看要将大总和列的哪些索引添加到您的子集中，以及是否值得在每个列中添加可比数据或不可配对数据案例.
function [accumout,accumleft]=getseparated(accum, bigcol, smallcol, m)
data=accum2mat(accum);
'indexing'
bigindex=bin2dec(num2str((data(:,bigcol)+1)/2))+1;
[~,bn]=size(bigcol);
[~,sn]=size(smallcol);
'Hamming weight'
b_ind=randperm(2^bn)'; %'
[~,I]=sort(sum((double(dec2bin(b_ind-1,bn))-48),2));
temp=zeros(2^bn,4+sn);
w=waitbar(0,'Processing');
for k=1:2^bn;
small_data=data(bigindex==b_ind(I(k)),smallcol);
if small_data
small_accum=mat2accum(small_data);
[small_accumpairs, small_accum]=getpairs(small_accum);
n_pairs=sum(small_accumpairs);
n_non_pairs=sum(small_accum);
sum_non_pairs=sum(accum2mat(small_accum));
else
n_pairs=0;
n_non_pairs=0;
sum_non_pairs=zeros(1,sn);
end
ham_weight=sum((double(dec2bin(b_ind(I(k))-1,bn))-48),2);
temp(k,:)=[b_ind(I(k)) n_pairs n_non_pairs ham_weight sum_non_pairs];
waitbar(k/2^bn);
end
close(w)
pair_ind=1;
nonpair_ind=1;
runningsum=[0 0 0 0 0 0 0 0 0 0];
temp2=zeros(2^bn,2);
while sum(sum(temp2))<=m
if pair_ind<=2^bn
pairsum=[(((double(dec2bin((temp(pair_ind,1)-1),bn))-48)*2)-1)*temp(pair_ind,2) zeros(1,sn)];
end
if nonpair_ind<=2^bn
nonpairsum=[(((double(dec2bin((temp(nonpair_ind,1)-1),bn))-48)*2)-1)*temp(nonpair_ind,3) temp(nonpair_ind,5:5+sn-1)];
end
if nonpair_ind==(2^bn)+1
temp2(pair_ind,1)=temp(pair_ind,2);
runningsum=runningsum+pairsum;
pair_ind=pair_ind+1;
elseif pair_ind==(2^bn)+1
temp2(nonpair_ind,2)=temp(nonpair_ind,3);
runningsum=runningsum+nonpairsum;
nonpair_ind=nonpair_ind+1;
elseif sum(abs(runningsum+pairsum))<=sum(abs(runningsum+nonpairsum))
temp2(pair_ind,1)=temp(pair_ind,2);
runningsum=runningsum+pairsum;
pair_ind=pair_ind+1;
elseif sum(abs(runningsum+pairsum))>sum(abs(runningsum+nonpairsum))
temp2(nonpair_ind,2)=temp(nonpair_ind,3);
runningsum=runningsum+nonpairsum;
nonpair_ind=nonpair_ind+1;
end
end
accumout=zeros(2^(bn+sn),1);
for k=1:2^bn
if temp2(k,:)
small_data=data(bigindex==temp(k,1),smallcol);
if small_data
small_accum=mat2accum(small_data);
[small_accumpairs, small_accum]=getpairs(small_accum);
pairs=accum2mat(small_accumpairs);
non_pairs=accum2mat(small_accum);
else
pairs=zeros(1,sn);
non_pairs=zeros(1,sn);
end
if temp2(k,1)
datatemp=zeros(temp2(k,1),sn+bn);
datatemp(:,bigcol)=((double(dec2bin(ones(temp2(k,1),1)*(temp(k,1)-1),bn))-48)*2)-1;
datatemp(:,smallcol)=pairs;
accumout=accumout+mat2accum(datatemp);
end
if temp2(k,2)
datatemp=zeros(temp2(k,2),sn+bn);
datatemp(:,bigcol)=((double(dec2bin(ones(temp2(k,2),1)*(temp(k,1)-1),bn))-48)*2)-1;
datatemp(:,smallcol)=non_pairs;
accumout=accumout+mat2accum(datatemp);
end
end
end
accumleft=accum-accumout;
对于由第一种情况的 5 列和第二种情况的 5 列组成的数据，可以构造一组 100k 行，小和列中的总和小于 1k，而小和列中的总和在 10k 到 30k 之间大的。
值得注意的是，数据的大小、所需子集的大小、1和-1的分布，都会对该算法的性能产生很大的影响ms。
由于您的回答似乎表明您有兴趣寻找更大的序列（更大的 n），下面的代码试图找到最大的 n，允许删除最多 10% 的行（即 25,000）。也就是说，代码通过从集合中删除最佳行最多 25,000 次来最小化完整数据集的 sum( abs( sum( data, 1)))。这应该与最小化均值（您陈述的问题）相同。该代码使用 [1, 1024] 范围内的索引来提高效率，直到最后一步产生最终输出。顺序变量设置为 10（您陈述的问题）对应于 2^10 = 1024 个可能的行向量。通过将所有 -1 值设置为 0 并采用二进制表示来找到给定行向量的索引，例如 [-1 -1 -1 -1 -1 -1 -1 -1 1]。所以在这个例子中，行向量的索引是 [0 0 0 0 0 0 0 0 0 1] = 1。 （请注意，索引 1 实际上已转换为 2，因为 MATLAB 不允许索引为 0。）
我已经测试了它的均匀随机分布（一个简单的案例），它通常在删除 ~1000 行后收敛到真实的最小值（即 sum( abs( sum( data, 1))) = 0）。 Click here to run the example code below for the uniform random case on AlgorithmHub。每次 运行 都会选择一个新的随机集，通常需要大约 30 秒才能完成该基础设施。
Click here to upload a csv file of your data set and run the example code on AlgorithmHub。 link 到 output.cvs 将允许您下载结果。如果您想获得特定的 n，应该很容易修改代码以支持您添加新行的方法。使用索引思想 should 和相应的查找 table (lut) 将有助于保持这种效率。否则，如果您想要一个特定的大 n，即使总和为 0（最小值），您也可以继续删除行。
% Generate data set as vector of length order with elements in set {1,-1}.
tic();
rows = 250000;
order = 10;
rowFraction = 0.1;
maxRowsToRemove = rows * rowFraction;
data = rand( rows, order);
data( data >= 0.5) = 1;
data( data < 0.5) = -1;
% Convert data to an index to one of 2^order vectors of 1 or -1.
% We set the -1 values to 0 and get the binary representation of the
% vector of binary values.
a = data;
a( a==-1)=0;
ndx = zeros(1,length(a));
ndx(:) = a(:,1)*2^9+a(:,2)*2^8+a(:,3)*2^7+a(:,4)*2^6+a(:,5)*2^5+...
a(:,6)*2^4+a(:,7)*2^3+a(:,8)*2^2+a(:,9)*2+a(:,10) + 1;
% Determine how many of each index we have in data pool.
bins = zeros( 1, 2^order);
binsRemoved = zeros( 1, 2^order);
for ii = 1:length( ndx)
bins( ndx(ii)) = bins( ndx(ii)) + 1;
end
colSum = sum(data,1);
sumOfColSum = sum(abs(colSum));
absSum = sumOfColSum;
lut = genLutForNdx( order);
nRemoved = 0;
curSum = colSum;
for ii = 1:maxRowsToRemove
if ( absSum == 0)
disp( sprintf( '\nminimum solution found'));
break;
end
ndxR = findNdxToRemove( curSum, bins, lut);
if ndxR > 0
bins( ndxR) = bins( ndxR) - 1;
binsRemoved( ndxR) = binsRemoved( ndxR) + 1;
curSum = curSum - lut( ndxR, :);
nRemoved = nRemoved + 1;
absSum = sum( abs( curSum));
else
disp( sprintf( '\nearly termination'));
break;
end
end
stat1 = sprintf( ...
'stats-L1: original sum = %d, final sum = %d, num rows removed = %d',...
sumOfColSum, absSum, nRemoved);
stat2 = sprintf( ...
'stats-L2: iter = %d, run time = %.2f sec\n', ii, toc());
disp( stat1);
disp( stat2);
% Show list of indicies removed along with the number of each removed.
binRndx = find( binsRemoved != 0);
ndxRemovedHist = [binRndx', binsRemoved(binRndx(:))'];
disp( sprintf( '%s\t%s', 'INDEX', 'NUM_REMOVED'));
for ii = 1: length( ndxRemovedHist)
disp( sprintf( '%d\t%d', ndxRemovedHist(ii,1), ndxRemovedHist(ii,2)));
end
% Generate the modified data array from the list of removed elements.
modData = data;
lr = [];
for ii = 1: length( ndxRemovedHist)
sr = find( ndx==ndxRemovedHist(ii,1));
lr = [lr, sr(1:ndxRemovedHist(ii,2))];
end
modData( lr, :) = [];
disp( sprintf( 'modified data array in variable "modData"'));
% ****************************************************
% Generate data set as vector of length order with elements in set {1,-1}.
tic();
rows = 250000;
order = 10;
rowFraction = 0.1;
maxRowsToRemove = rows * rowFraction;
data = rand( rows, order);
data( data >= 0.5) = 1;
data( data < 0.5) = -1;
% Convert data to an index to one of 2^order vectors of 1 or -1.
% We set the -1 values to 0 and get the binary representation of the
% vector of binary values.
a = data;
a( a==-1)=0;
ndx = zeros(1,length(a));
ndx(:) = a(:,1)*2^9+a(:,2)*2^8+a(:,3)*2^7+a(:,4)*2^6+a(:,5)*2^5+...
a(:,6)*2^4+a(:,7)*2^3+a(:,8)*2^2+a(:,9)*2+a(:,10) + 1;
% Determine how many of each index we have in data pool.
bins = zeros( 1, 2^order);
binsRemoved = zeros( 1, 2^order);
for ii = 1:length( ndx)
bins( ndx(ii)) = bins( ndx(ii)) + 1;
end
colSum = sum(data,1);
sumOfColSum = sum(abs(colSum));
absSum = sumOfColSum;
lut = genLutForNdx( order);
nRemoved = 0;
curSum = colSum;
for ii = 1:maxRowsToRemove
if ( absSum == 0)
disp( sprintf( '\nminimum solution found'));
break;
end
ndxR = findNdxToRemove( curSum, bins, lut);
if ndxR > 0
bins( ndxR) = bins( ndxR) - 1;
binsRemoved( ndxR) = binsRemoved( ndxR) + 1;
curSum = curSum - lut( ndxR, :);
nRemoved = nRemoved + 1;
absSum = sum( abs( curSum));
else
disp( sprintf( '\nearly termination'));
break;
end
end
stat1 = sprintf( ...
'stats-L1: original sum = %d, final sum = %d, num rows removed = %d',...
sumOfColSum, absSum, nRemoved);
stat2 = sprintf( ...
'stats-L2: iter = %d, run time = %.2f sec\n', ii, toc());
disp( stat1);
disp( stat2);
% Show list of indicies removed along with the number of each removed.
binRndx = find( binsRemoved != 0);
ndxRemovedHist = [binRndx', binsRemoved(binRndx(:))'];
disp( sprintf( '%s\t%s', 'INDEX', 'NUM_REMOVED'));
for ii = 1: length( ndxRemovedHist)
disp( sprintf( '%d\t%d', ndxRemovedHist(ii,1), ndxRemovedHist(ii,2)));
end
% Generate the modified data array from the list of removed elements.
modData = data;
lr = [];
for ii = 1: length( ndxRemovedHist)
sr = find( ndx==ndxRemovedHist(ii,1));
lr = [lr, sr(1:ndxRemovedHist(ii,2))];
end
modData( lr, :) = [];
disp( sprintf( 'modified data array in variable "modData"'));
% ****************************************************
function ndx = findNdxToRemove( curSum, bins, lut)
% See if ideal index to remove exists in current bin set. We look at the
% sign of each element of the current sum to determine index to remove
aa = zeros( size( curSum));
if (isempty( find( curSum == 0)))
aa( curSum < 0) = 0;
aa( curSum > 0) = 1;
ndx = aa(1)*2^9+aa(2)*2^8+aa(3)*2^7+aa(4)*2^6+aa(5)*2^5+...
aa(6)*2^4+aa(7)*2^3+aa(8)*2^2+aa(9)*2+aa(10) + 1;
if( bins(ndx) > 0)
% Optimal row to remove was found directly.
return;
end
end
% Serach through all the non-empty indices that remain for best to remove.
delta = 0;
ndx = -1;
minSum = sum( abs( curSum));
minSumOrig = minSum;
bestNdx = -1;
firstFound = 1;
for ii = 1:length( bins)
if ( bins(ii) > 0)
tmp = sum( abs( curSum - lut( ii,:)));
if ( firstFound)
minSum = tmp;
bestNdx = ii;
firstFound = 0;
elseif ( tmp < minSum)
minSum = tmp;
bestNdx = ii;
end
end
end
ndx = bestNdx;
下面的方法呢。由于 10 列只有 +1 和 -1 值，所以只有 1024 行可能。所以我们现在的数据是：
- 具有 -1 和 +1 系数的 1024 x 10 矩阵
a(i,j)。该矩阵具有所有不同的可能（唯一）行。
- 一个向量
v(i)，其中包含我们看到第 i 行的次数。
现在我们可以写一个简单的混合整数规划问题如下：
备注：
- 我们只有1024个整型变量
- 我们在 x(i) 上设置了一个上限，表示可以选择一行的次数
- 我们使用所谓的变量拆分技术对绝对值建模并保持模型线性
- 最小化均值和最小化和是一样的（区别是常数因子）
- 关于 optcr 的行告诉 MIP 求解器找到经过验证的全局最优解
- 一个好的 MIP 求解器应该能够很快找到解决方案。我使用 250k 行和 N=100 对一些随机数据进行了测试。
我实际上认为这是一个简单的问题。
- 重申一下：此方法提供了经过验证的全局最优解。
- 可以找到更多详细信息 here。
我有一个大数组（大约 250,000 x 10）。每行包含 1 或 -1。例如：
data(1, :) = [1, -1, -1, -1, -1, -1, -1, -1, 1, -1];
我需要 select 组 n 行，这样 列的绝对总和的平均值 被最小化（尽可能接近零）。因此，在这个玩具示例中，其中 n=2：
[ 1 1 1 1]
[-1 -1 -1 -1]
[-1 1 -1 1]
我会 select 第 1 行和第 2 行，因为它们总和为 [0 0 0 0]（平均值 0），这是 n=2 时可能的最小值。
我尝试了下面建议的方法（寻找互补对），但对于我的数据集，这只能形成 23k 行的平衡子集。因此，我需要一个近似值，它生成大小为 n 行的子集，但具有列的绝对和的最小均值。
到目前为止我发现的最佳方法如下：选择一个起始子集，将剩余的每一行迭代添加到基数，如果它提高了列的绝对和的平均值则保留它。这是非常粗糙的，我相信有更好的方法。它也很容易卡在错误的最小值，所以需要添加一个意外事件：
shuffle = randperm(size(data));
data_shuffled = data(shuffle, :);
base = data_shuffled(1:30000, :);
pool = data_shuffled(30001:end, :);
best_mean = mean(abs(sum(base, 1)));
best_matrix = base;
n = 100000;
for k = 1:20
for i = 1:size(pool, 1)
temp = pool (i, :);
if(~isnan(temp(1)))
temp_sum = sum(temp, 1);
new_sum = temp_sum + sum(best, 1);
temp_mean = mean(abs(new_sum));
if(temp_mean < best_mean)
best_mean = temp_mean;
best_matrix = vertcat(best_matrix, temp);
pool(i, :) = NaN(1, 10);
end
end
end
if(size(best_matrix, 1) > n)
return
end
end
这实现了 ~17,000 列的绝对总和的平均值，这还不错。重复使用不同的种子可能会有所改善。
理想情况下，我不会将新元素添加到 best_matrix 的末尾，而是 将它 替换为一些元素，以实现最佳改进。
更新：我不想给出数据集的具体细节，因为所有解决方案都应该适用于指定格式的任何矩阵。
感谢大家的贡献！
不幸的是，这个问题超出了常规（连续）优化的范围。你的问题，可以参数化如下：
min_{S∈S_n} Σ_{j∈S}|Σ_i data_ji|
其中S_n是索引j∈{0,...,250000}的n个元素组合的集合，也可以重写为非常类似x中的正则二次整数规划问题：
min_x x'* data *data' *x
0<=x<=1 and x*1=n
其中 data 是您的 250000*10 矩阵，x 是我们正在寻找的 250000*1 组合向量。 （现在我们优化的是平方和，而不是绝对值的和...）
这道题被证明是NP-hard，也就是说要求全局最小值，必须遍历n种所有可能的组合，抽取250000种可能，这等于二项式系数（250000 n），等于250000!/(n!*(250000-n)!)...
祝你好运...;)
编辑
如果您要启发式地解决这个问题，因为我想您将需要一个解决方案，请使用启发式 here 而不是您的方法。
正如其他人所说，最佳解决方案可能是不可能的，所以我将专注于特定情况。
首先我假设每列的分布是独立的。
然后我在累加器 space 上工作，以减少数据大小并加快代码速度。
我通过将每个 -1 视为 0 并将每一行视为一个数字，并添加 1 以避免使用 0 作为索引，例如：
data(1,:)=[-1 1 -1 1 -1 1 -1 1 -1 1] -> '0101010101' -> 341 -> 342
有了这个我们可以累积数据为：
function accum=mat2accum(data)
[~,n]=size(data);
indexes=bin2dec(num2str((data+1)/2))+1;
accum=accumarray(indexes,1,[2^n 1]);
我考虑的第一种情况是当每列的总和与数据的大小相比较小时，这意味着所有列中的 1 和 -1 的数量相似。
sum(data) << size(data)
对于这种情况，您可以找到所有相互抵消的对，例如：
data(1,:)=[-1 1 -1 1 -1 1 -1 1 -1 1] -> '0101010101' -> 341 -> 342
data(2,:)=[1 -1 1 -1 1 -1 1 -1 1 -1] -> '1010101010' -> 682 -> 683
并且我们知道每一对都将位于累加器索引中的镜像位置，因此我们可以通过以下方式获得所有可能的对：
function [accumpairs, accumleft]=getpairs(accum)
accumpairs=min([accum,accum(end:-1:1)],[],2);
accumleft=accum-accumpairs;
使用随机生成的数据，我能够在一组 250k 行中获得 >100k 对，并且每对的子集在每列中的总和等于零。因此，如果 1 和 -1 均匀分布，这可能就足够了。
我考虑的第二种情况是每列的总和远不为零，这意味着 1 和 -1 之间存在很大的不成比例。
abs(sum(data)) >> 0
通过反转总和为负的每一列，这不会影响数据，因为最后可以再次反转这些列。可以强制不成比例为 1 多于 -1。并且通过首先提取该数据的可能对，不成比例更加明显。
有了这样准备的数据，就可以将问题视为最小化所需集合中 1 的数量。为此，我们首先将可能的指标随机化，然后计算并排序每个指标的汉明权重（二进制表示中1的个数），然后收集汉明权重尽可能小的数据。
function [accumlast,accumleft]=resto(accum,m)
[N,~]=size(accum);
columns=log2(N);
indexes=randperm(N)'; %'
[~,I]=sort(sum((double(dec2bin(indexes-1,columns))-48),2));
accumlast=zeros(N,1);
for k=indexes(I)' %'
accumlast(k)=accum(k);
if sum(accumlast)>=m
break
end
end
accumleft=accum-accumlast;
使用随机生成的数据，其中 1 的数量大约是 -1 的 2 倍，并且每列的总和约为 80k，我可以找到 100k 数据的子集，每列的总和约为 5k。
第三种情况，有些列的总和接近于零，有些则不是。在这种情况下，您将列分为具有大和的列和具有小和的列，然后按大和列的汉明权重对数据进行排序，并在每个大列索引中获取小和列对.这将为大和列的每个索引创建一个矩阵，其中包含可能对的数量、不成对的行数以及小列的不成对行的总和。
现在您可以使用该信息来保存 运行 总和，并查看要将大总和列的哪些索引添加到您的子集中，以及是否值得在每个列中添加可比数据或不可配对数据案例.
function [accumout,accumleft]=getseparated(accum, bigcol, smallcol, m)
data=accum2mat(accum);
'indexing'
bigindex=bin2dec(num2str((data(:,bigcol)+1)/2))+1;
[~,bn]=size(bigcol);
[~,sn]=size(smallcol);
'Hamming weight'
b_ind=randperm(2^bn)'; %'
[~,I]=sort(sum((double(dec2bin(b_ind-1,bn))-48),2));
temp=zeros(2^bn,4+sn);
w=waitbar(0,'Processing');
for k=1:2^bn;
small_data=data(bigindex==b_ind(I(k)),smallcol);
if small_data
small_accum=mat2accum(small_data);
[small_accumpairs, small_accum]=getpairs(small_accum);
n_pairs=sum(small_accumpairs);
n_non_pairs=sum(small_accum);
sum_non_pairs=sum(accum2mat(small_accum));
else
n_pairs=0;
n_non_pairs=0;
sum_non_pairs=zeros(1,sn);
end
ham_weight=sum((double(dec2bin(b_ind(I(k))-1,bn))-48),2);
temp(k,:)=[b_ind(I(k)) n_pairs n_non_pairs ham_weight sum_non_pairs];
waitbar(k/2^bn);
end
close(w)
pair_ind=1;
nonpair_ind=1;
runningsum=[0 0 0 0 0 0 0 0 0 0];
temp2=zeros(2^bn,2);
while sum(sum(temp2))<=m
if pair_ind<=2^bn
pairsum=[(((double(dec2bin((temp(pair_ind,1)-1),bn))-48)*2)-1)*temp(pair_ind,2) zeros(1,sn)];
end
if nonpair_ind<=2^bn
nonpairsum=[(((double(dec2bin((temp(nonpair_ind,1)-1),bn))-48)*2)-1)*temp(nonpair_ind,3) temp(nonpair_ind,5:5+sn-1)];
end
if nonpair_ind==(2^bn)+1
temp2(pair_ind,1)=temp(pair_ind,2);
runningsum=runningsum+pairsum;
pair_ind=pair_ind+1;
elseif pair_ind==(2^bn)+1
temp2(nonpair_ind,2)=temp(nonpair_ind,3);
runningsum=runningsum+nonpairsum;
nonpair_ind=nonpair_ind+1;
elseif sum(abs(runningsum+pairsum))<=sum(abs(runningsum+nonpairsum))
temp2(pair_ind,1)=temp(pair_ind,2);
runningsum=runningsum+pairsum;
pair_ind=pair_ind+1;
elseif sum(abs(runningsum+pairsum))>sum(abs(runningsum+nonpairsum))
temp2(nonpair_ind,2)=temp(nonpair_ind,3);
runningsum=runningsum+nonpairsum;
nonpair_ind=nonpair_ind+1;
end
end
accumout=zeros(2^(bn+sn),1);
for k=1:2^bn
if temp2(k,:)
small_data=data(bigindex==temp(k,1),smallcol);
if small_data
small_accum=mat2accum(small_data);
[small_accumpairs, small_accum]=getpairs(small_accum);
pairs=accum2mat(small_accumpairs);
non_pairs=accum2mat(small_accum);
else
pairs=zeros(1,sn);
non_pairs=zeros(1,sn);
end
if temp2(k,1)
datatemp=zeros(temp2(k,1),sn+bn);
datatemp(:,bigcol)=((double(dec2bin(ones(temp2(k,1),1)*(temp(k,1)-1),bn))-48)*2)-1;
datatemp(:,smallcol)=pairs;
accumout=accumout+mat2accum(datatemp);
end
if temp2(k,2)
datatemp=zeros(temp2(k,2),sn+bn);
datatemp(:,bigcol)=((double(dec2bin(ones(temp2(k,2),1)*(temp(k,1)-1),bn))-48)*2)-1;
datatemp(:,smallcol)=non_pairs;
accumout=accumout+mat2accum(datatemp);
end
end
end
accumleft=accum-accumout;
对于由第一种情况的 5 列和第二种情况的 5 列组成的数据，可以构造一组 100k 行，小和列中的总和小于 1k，而小和列中的总和在 10k 到 30k 之间大的。
值得注意的是，数据的大小、所需子集的大小、1和-1的分布，都会对该算法的性能产生很大的影响ms。
由于您的回答似乎表明您有兴趣寻找更大的序列（更大的 n），下面的代码试图找到最大的 n，允许删除最多 10% 的行（即 25,000）。也就是说，代码通过从集合中删除最佳行最多 25,000 次来最小化完整数据集的 sum( abs( sum( data, 1)))。这应该与最小化均值（您陈述的问题）相同。该代码使用 [1, 1024] 范围内的索引来提高效率，直到最后一步产生最终输出。顺序变量设置为 10（您陈述的问题）对应于 2^10 = 1024 个可能的行向量。通过将所有 -1 值设置为 0 并采用二进制表示来找到给定行向量的索引，例如 [-1 -1 -1 -1 -1 -1 -1 -1 1]。所以在这个例子中，行向量的索引是 [0 0 0 0 0 0 0 0 0 1] = 1。 （请注意，索引 1 实际上已转换为 2，因为 MATLAB 不允许索引为 0。）
我已经测试了它的均匀随机分布（一个简单的案例），它通常在删除 ~1000 行后收敛到真实的最小值（即 sum( abs( sum( data, 1))) = 0）。 Click here to run the example code below for the uniform random case on AlgorithmHub。每次 运行 都会选择一个新的随机集，通常需要大约 30 秒才能完成该基础设施。
Click here to upload a csv file of your data set and run the example code on AlgorithmHub。 link 到 output.cvs 将允许您下载结果。如果您想获得特定的 n，应该很容易修改代码以支持您添加新行的方法。使用索引思想 should 和相应的查找 table (lut) 将有助于保持这种效率。否则，如果您想要一个特定的大 n，即使总和为 0（最小值），您也可以继续删除行。
% Generate data set as vector of length order with elements in set {1,-1}.
tic();
rows = 250000;
order = 10;
rowFraction = 0.1;
maxRowsToRemove = rows * rowFraction;
data = rand( rows, order);
data( data >= 0.5) = 1;
data( data < 0.5) = -1;
% Convert data to an index to one of 2^order vectors of 1 or -1.
% We set the -1 values to 0 and get the binary representation of the
% vector of binary values.
a = data;
a( a==-1)=0;
ndx = zeros(1,length(a));
ndx(:) = a(:,1)*2^9+a(:,2)*2^8+a(:,3)*2^7+a(:,4)*2^6+a(:,5)*2^5+...
a(:,6)*2^4+a(:,7)*2^3+a(:,8)*2^2+a(:,9)*2+a(:,10) + 1;
% Determine how many of each index we have in data pool.
bins = zeros( 1, 2^order);
binsRemoved = zeros( 1, 2^order);
for ii = 1:length( ndx)
bins( ndx(ii)) = bins( ndx(ii)) + 1;
end
colSum = sum(data,1);
sumOfColSum = sum(abs(colSum));
absSum = sumOfColSum;
lut = genLutForNdx( order);
nRemoved = 0;
curSum = colSum;
for ii = 1:maxRowsToRemove
if ( absSum == 0)
disp( sprintf( '\nminimum solution found'));
break;
end
ndxR = findNdxToRemove( curSum, bins, lut);
if ndxR > 0
bins( ndxR) = bins( ndxR) - 1;
binsRemoved( ndxR) = binsRemoved( ndxR) + 1;
curSum = curSum - lut( ndxR, :);
nRemoved = nRemoved + 1;
absSum = sum( abs( curSum));
else
disp( sprintf( '\nearly termination'));
break;
end
end
stat1 = sprintf( ...
'stats-L1: original sum = %d, final sum = %d, num rows removed = %d',...
sumOfColSum, absSum, nRemoved);
stat2 = sprintf( ...
'stats-L2: iter = %d, run time = %.2f sec\n', ii, toc());
disp( stat1);
disp( stat2);
% Show list of indicies removed along with the number of each removed.
binRndx = find( binsRemoved != 0);
ndxRemovedHist = [binRndx', binsRemoved(binRndx(:))'];
disp( sprintf( '%s\t%s', 'INDEX', 'NUM_REMOVED'));
for ii = 1: length( ndxRemovedHist)
disp( sprintf( '%d\t%d', ndxRemovedHist(ii,1), ndxRemovedHist(ii,2)));
end
% Generate the modified data array from the list of removed elements.
modData = data;
lr = [];
for ii = 1: length( ndxRemovedHist)
sr = find( ndx==ndxRemovedHist(ii,1));
lr = [lr, sr(1:ndxRemovedHist(ii,2))];
end
modData( lr, :) = [];
disp( sprintf( 'modified data array in variable "modData"'));
% ****************************************************
% Generate data set as vector of length order with elements in set {1,-1}.
tic();
rows = 250000;
order = 10;
rowFraction = 0.1;
maxRowsToRemove = rows * rowFraction;
data = rand( rows, order);
data( data >= 0.5) = 1;
data( data < 0.5) = -1;
% Convert data to an index to one of 2^order vectors of 1 or -1.
% We set the -1 values to 0 and get the binary representation of the
% vector of binary values.
a = data;
a( a==-1)=0;
ndx = zeros(1,length(a));
ndx(:) = a(:,1)*2^9+a(:,2)*2^8+a(:,3)*2^7+a(:,4)*2^6+a(:,5)*2^5+...
a(:,6)*2^4+a(:,7)*2^3+a(:,8)*2^2+a(:,9)*2+a(:,10) + 1;
% Determine how many of each index we have in data pool.
bins = zeros( 1, 2^order);
binsRemoved = zeros( 1, 2^order);
for ii = 1:length( ndx)
bins( ndx(ii)) = bins( ndx(ii)) + 1;
end
colSum = sum(data,1);
sumOfColSum = sum(abs(colSum));
absSum = sumOfColSum;
lut = genLutForNdx( order);
nRemoved = 0;
curSum = colSum;
for ii = 1:maxRowsToRemove
if ( absSum == 0)
disp( sprintf( '\nminimum solution found'));
break;
end
ndxR = findNdxToRemove( curSum, bins, lut);
if ndxR > 0
bins( ndxR) = bins( ndxR) - 1;
binsRemoved( ndxR) = binsRemoved( ndxR) + 1;
curSum = curSum - lut( ndxR, :);
nRemoved = nRemoved + 1;
absSum = sum( abs( curSum));
else
disp( sprintf( '\nearly termination'));
break;
end
end
stat1 = sprintf( ...
'stats-L1: original sum = %d, final sum = %d, num rows removed = %d',...
sumOfColSum, absSum, nRemoved);
stat2 = sprintf( ...
'stats-L2: iter = %d, run time = %.2f sec\n', ii, toc());
disp( stat1);
disp( stat2);
% Show list of indicies removed along with the number of each removed.
binRndx = find( binsRemoved != 0);
ndxRemovedHist = [binRndx', binsRemoved(binRndx(:))'];
disp( sprintf( '%s\t%s', 'INDEX', 'NUM_REMOVED'));
for ii = 1: length( ndxRemovedHist)
disp( sprintf( '%d\t%d', ndxRemovedHist(ii,1), ndxRemovedHist(ii,2)));
end
% Generate the modified data array from the list of removed elements.
modData = data;
lr = [];
for ii = 1: length( ndxRemovedHist)
sr = find( ndx==ndxRemovedHist(ii,1));
lr = [lr, sr(1:ndxRemovedHist(ii,2))];
end
modData( lr, :) = [];
disp( sprintf( 'modified data array in variable "modData"'));
% ****************************************************
function ndx = findNdxToRemove( curSum, bins, lut)
% See if ideal index to remove exists in current bin set. We look at the
% sign of each element of the current sum to determine index to remove
aa = zeros( size( curSum));
if (isempty( find( curSum == 0)))
aa( curSum < 0) = 0;
aa( curSum > 0) = 1;
ndx = aa(1)*2^9+aa(2)*2^8+aa(3)*2^7+aa(4)*2^6+aa(5)*2^5+...
aa(6)*2^4+aa(7)*2^3+aa(8)*2^2+aa(9)*2+aa(10) + 1;
if( bins(ndx) > 0)
% Optimal row to remove was found directly.
return;
end
end
% Serach through all the non-empty indices that remain for best to remove.
delta = 0;
ndx = -1;
minSum = sum( abs( curSum));
minSumOrig = minSum;
bestNdx = -1;
firstFound = 1;
for ii = 1:length( bins)
if ( bins(ii) > 0)
tmp = sum( abs( curSum - lut( ii,:)));
if ( firstFound)
minSum = tmp;
bestNdx = ii;
firstFound = 0;
elseif ( tmp < minSum)
minSum = tmp;
bestNdx = ii;
end
end
end
ndx = bestNdx;
下面的方法呢。由于 10 列只有 +1 和 -1 值，所以只有 1024 行可能。所以我们现在的数据是：
- 具有 -1 和 +1 系数的 1024 x 10 矩阵
a(i,j)。该矩阵具有所有不同的可能（唯一）行。 - 一个向量
v(i)，其中包含我们看到第 i 行的次数。
现在我们可以写一个简单的混合整数规划问题如下：
备注：
- 我们只有1024个整型变量
- 我们在 x(i) 上设置了一个上限，表示可以选择一行的次数
- 我们使用所谓的变量拆分技术对绝对值建模并保持模型线性
- 最小化均值和最小化和是一样的（区别是常数因子）
- 关于 optcr 的行告诉 MIP 求解器找到经过验证的全局最优解
- 一个好的 MIP 求解器应该能够很快找到解决方案。我使用 250k 行和 N=100 对一些随机数据进行了测试。 我实际上认为这是一个简单的问题。
- 重申一下：此方法提供了经过验证的全局最优解。
- 可以找到更多详细信息 here。