arrays - una - Minimizando la suma de columnas de matriz en Matlab

seleccionar elementos de una matriz matlab

Question

Tengo una gran variedad (aprox. 250,000 x 10). Cada fila contiene 1s o -1s. P.ej:

data(1, :) = [1, -1, -1, -1, -1, -1, -1, -1, 1, -1];

Necesito seleccionar conjuntos de n filas, de modo que la media de las sumas absolutas de las columnas se minimice (lo más cerca posible de cero). Entonces, en este ejemplo de juguete, donde n = 2:

[ 1 1 1 1]
[-1 -1 -1 -1]
[-1 1 -1 1]

Seleccionaría las filas 1 y 2, ya que suman a [0 0 0 0] (media 0), que es el mínimo posible cuando n = 2.

Probé el método sugerido a continuación (para encontrar pares complementarios), pero para mi conjunto de datos, esto solo puede formar un subconjunto equilibrado de 23k filas. Por lo tanto, necesito una aproximación que genere un subconjunto de filas de tamaño n, pero con los medios mínimos de sumas absolutas de las columnas.

El mejor enfoque que he encontrado hasta ahora es el siguiente: elija un subconjunto de inicio, agregue de forma iterativa cada fila del resto a la base y manténgala si mejora la media de las sumas absolutas de las columnas. Esto es muy crudo y estoy seguro de que hay mejores maneras. También es propenso a atascarse en falsos mínimos, por lo que se debe agregar una contingencia:

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

Esto logra una media de las sumas absolutas de las columnas de ~ 17,000, lo cual no es tan malo. Repetir con diferentes semillas probablemente lo mejore un poco.

Idealmente, en lugar de simplemente agregar el nuevo elemento al final de best_matrix, lo cambiaría por algún elemento, a fin de lograr la mejor mejora.

Actualización: no quiero dar detalles específicos del conjunto de datos porque todas las soluciones deben ser aplicables a cualquier matriz en el formato especificado.

¡Gracias a todos los que contribuyeron!

Answer 1

¿Qué pasa con el siguiente enfoque. Con 10 columnas que solo tienen valores de +1 y -1, solo son posibles 1024 filas diferentes. Así que nuestros datos son ahora:

1. a 1024 x 10 matriz a(i,j) con -1 y +1 coeficientes. Esta matriz tiene todas las diferentes filas posibles (únicas).
2. un vector v(i) con el número de veces que vimos la fila i.

Ahora podemos escribir un problema simple de programación de enteros mixtos de la siguiente manera:

Notas:

• Solo tenemos 1024 variables enteras.
• Establecemos un límite superior en x (i) que indica cuántas veces se puede seleccionar una fila
• Utilizamos una técnica llamada de división de variables para modelar los valores absolutos y mantener el modelo lineal.
• Minimizar la media es lo mismo que minimizar la suma (la diferencia es un factor constante)
• La línea sobre optcr le dice al solucionador de MIP que encuentre soluciones óptimas globales probadas
• Un buen solucionador de MIP debería poder encontrar soluciones muy rápidamente. Probé con algunos datos aleatorios usando 250k filas y N = 100. Realmente creo que este es un problema fácil.
• Para reiterar: este método ofrece soluciones óptimas globales comprobadas.
• Algunos detalles más se pueden encontrar here .

Answer 2

Como otros habían dicho que una solución óptima podría ser imposible, me centraré en casos específicos.

Primero asumo independencia de las distribuciones de cada columna.

Luego trabajo en el espacio del acumulador para reducir el tamaño de los datos y acelerar el código.

Lo hago tomando cada -1 como 0 y considerando cada fila como un número, y agrego 1 para evitar trabajar con 0 como un índice, como:

data(1,:)=[-1 1 -1 1 -1 1 -1 1 -1 1] -> ''0101010101'' -> 341 -> 342

Con esto podemos acumular los datos como:

function accum=mat2accum(data)
[~,n]=size(data);
indexes=bin2dec(num2str((data+1)/2))+1;
accum=accumarray(indexes,1,[2^n 1]);

El primer caso que considero es cuando la suma de cada columna es un número pequeño en comparación con el tamaño de los datos, esto significa que hay una cantidad similar de 1 y -1 en todas las columnas.

sum(data) << size(data)

Para este caso puedes encontrar todas las parejas que se cancelan entre sí como:

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

Y sabemos que cada par estará en la posición reflejada en el índice del acumulador, por lo que podemos obtener todos los pares posibles con:

function [accumpairs, accumleft]=getpairs(accum)
accumpairs=min([accum,accum(end:-1:1)],[],2);
accumleft=accum-accumpairs;

Con los datos generados al azar pude obtener> 100k pares en un conjunto de 250k filas, y un subconjunto de pares tendrá una suma igual a cero en cada columna. Entonces, si 1 y -1 están distribuidos equitativamente, esto podría ser suficiente.

El segundo caso que consideré fue cuando la suma de cada columna está lejos de cero, lo que significa que hay una gran desproporción entre 1 y -1.

abs(sum(data)) >> 0

Al invertir cada columna donde la suma es negativa, esto no afectaría a los datos, ya que al final es posible invertir esas columnas nuevamente. Es posible forzar que la desproporción sea más de 1 que de -1. Y al extraer primero los posibles pares de estos datos, la desproporción es aún más pronunciada.

Con los datos preparados como tales, es posible tratar el problema para minimizar el número de 1 en el conjunto requerido. Para esto, primero aleatorizamos los índices posibles, luego calculamos y clasificamos el peso de Hamming (número de 1 en la representación binaria) de cada índice, y luego recopilamos los datos con el peso de Hamming más pequeño posible.

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;

Con datos generados aleatoriamente donde hubo aproximadamente 2 veces más 1 que -1, y la suma de cada columna fue alrededor de 80k, puedo encontrar un subconjunto de datos de 100k con una suma de alrededor de 5k en cada columna.

El tercer caso, es cuando algunas columnas suman cerca de cero y otras no. En este caso, se separan las columnas en las que tienen una gran suma y las que tienen una pequeña suma, luego se clasifican los datos por el peso de Hamming de las columnas de la gran suma y se obtienen los pares de las columnas de la pequeña suma dentro de cada una de las columnas grandes. . Esto creará una matriz con el número de pares posibles, el número de filas no emparejadas y la suma de las filas no emparejadas de las columnas pequeñas, para cada índice de las columnas de suma grande.

Ahora puede usar esa información para mantener una suma corriente y ver qué índices de las columnas de la gran suma agregar a su subconjunto, y también si vale la pena agregar la parábola o los datos no pareables en cada caso.

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;

Con los datos compuestos por 5 columnas del primer caso y 5 columnas del segundo caso, fue posible construir un conjunto de 100k filas con <1k de suma en las columnas de suma pequeña y entre 10k y 30k en las grandes.

Vale la pena señalar que el tamaño de los datos, el tamaño del subconjunto requerido y la distribución de 1 y -1 tendrán un gran efecto en el rendimiento de estos algoritmos.

Answer 3

Lamentablemente, este problema queda fuera del alcance de la optimización regular (continua). Su problema, que se puede parametrizar de la siguiente manera:

min_{S∈S_n} Σ_{j∈S}|Σ_i data_ji|

Donde S_n es el conjunto de n-elementos combinaciones de índices j∈{0,...,250000} , también se puede reescribir como un problema de programación de enteros cuadráticos regular muy similar en x :

min_x x''* data *data'' *x
0<=x<=1 and x*1=n

Donde los data son su matriz de 250000 * 10 y x es el vector de combinaciones de 250000 * 1 que estamos buscando. (Ahora optimizamos la suma de cuadrados en lugar de la suma de valores absolutos ...)

Se demostró que este problema es NP-hard , lo que significa que para encontrar el minimizador global, debe pasar por todas las combinaciones posibles de n dibujados en 250000 posibilidades, que es igual al coeficiente binomial (250000 n), ¡que es igual a 250000!/(n!*(250000-n)!) ...

Buena suerte... ;)

EDITAR

Si va a resolver esto heurísticamente, ya que supongo que va a necesitar una solución, utilice la heurística here lugar de su enfoque.

Answer 4

Ya que sus respuestas parecían indicar que estaba interesado en encontrar secuencias más grandes (n más grande), el código siguiente intenta encontrar la n más grande que permita eliminar hasta el 10% de las filas (es decir, 25,000). Ese es el código que minimiza la sum( abs( sum( data, 1))) del conjunto completo de datos al eliminar la mejor fila del conjunto hasta 25,000 veces. Esto debería ser lo mismo que minimizar la media (su problema declarado). El código utiliza índices en el rango [1, 1024] por eficiencia hasta que se produce la salida final en el último paso. La variable de orden se establece en 10 (su problema declarado) correspondiente a 2^10 = 1024 vectores de fila posibles. Se encuentra un índice para un vector de fila dado, por ejemplo [-1 -1 -1 -1 -1 -1 -1 -1 1] , configurando todos los valores -1 a 0 y tomando la representación binaria. Entonces, en este ejemplo, el índice del vector de la fila es [0 0 0 0 0 0 0 0 0 1] = 1 . (Tenga en cuenta que el índice 1 se convierte realmente a 2 ya que MATLAB no permite un índice de 0).

He probado esto para una distribución aleatoria uniforme (un caso fácil) y generalmente converge a un mínimo verdadero (es decir, sum( abs( sum( data, 1))) = 0 ) después de eliminar ~ 1000 filas. Haga clic aquí para ejecutar el código de ejemplo a continuación para el caso aleatorio uniforme en AlgorithmHub . Se elegirá un nuevo conjunto aleatorio cada vez que se ejecute y, en general, tomará unos 30 segundos completar en esa infraestructura.

Haga clic aquí para cargar un archivo csv de su conjunto de datos y ejecute el código de ejemplo en AlgorithmHub . Un enlace a output.cvs te permitirá descargar los resultados. El código debe modificarse fácilmente para admitir su método de agregar nuevas filas si desea obtener un n específico. El uso de la idea de índice con la tabla de consulta correspondiente (lut) ayudará a mantener esto eficiente. De lo contrario, si desea una n grande específica, puede seguir eliminando filas incluso cuando la suma sea 0 (mínimo).

% 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;