matlab - una - Encontrar columnas proporcionales en matriz

seleccionar elementos de una matriz matlab

Question

Tengo una gran matriz (1,000 filas y 50,000 columnas). Sé que algunas columnas están correlacionadas (el rango es solo 100) y sospecho que algunas columnas son incluso proporcionales. ¿Cómo puedo encontrar tales columnas proporcionales? (Una forma sería ejecutar corr(M(:,j),M(:,k)) ), pero ¿hay algo más eficiente?

Answer 1

Si entiendo su problema correctamente, desea determinar las columnas en su matriz que son linealmente dependientes , lo que significa que una columna es proporcional o escalar múltiplo de otra. Hay un algoritmo muy básico basado en la Descomposición QR . Para la descomposición QR, puede tomar cualquier matriz y descomponerla en un producto de dos matrices: Q y R En otras palabras:

A = Q*R

Q es una matriz ortogonal con cada columna como un vector unitario, de modo que al multiplicar Q por su transposición se obtiene la matriz de identidad ( Q^{T}*Q = I ). R es una matriz triangular derecha o triangular superior. Una teoría muy útil de Golub y Van Loan en su libro de 1996: Matriz de Cómputos es que una matriz se considera de rango completo si todos los valores de los elementos diagonales de R no son cero. Debido a la precisión del punto flotante en las computadoras, tendremos que establecer un umbral y verificar si hay valores en la diagonal de R que sean mayores que esta tolerancia. Si lo es, entonces esta columna correspondiente es una columna independiente. Simplemente podemos encontrar el valor absoluto de todas las diagonales, luego verificamos si son mayores que algunas tolerancias.

Podemos modificar esto ligeramente para que busquemos valores que sean menores que la tolerancia, lo que significaría que la columna no es independiente . La forma en que llamaría la factorización QR es de esta manera:

[Q,R] = qr(A, 0);

Q y R son de lo que acabo de hablar, y usted especifica la matriz A como entrada. El segundo parámetro 0 significa producir una versión económica de Q y R , donde si esta matriz fuera rectangular (como en su caso), esto devolvería una matriz cuadrada donde las dimensiones son las más grandes de los dos tamaños. En otras palabras, si tuviera una matriz como 5 x 8 , producir una matriz de tamaño económico le dará una salida de 5 x 8 , mientras que si no especifica el 0 , obtendrá una matriz de 8 x 8 .

Ahora, lo que realmente necesitamos es este estilo de invocación:

[Q,R,E] = qr(A, 0);

En este caso, E sería un vector de permutación, tal que:

A(:,E) = Q*R;

La razón por la cual esto es útil es porque ordena las columnas de Q y R de tal manera que la primera columna de la versión reordenada es la más probable que es independiente, seguida por esas columnas en orden decreciente de "fuerza" . Como tal, E le diría cuán probable es que cada columna sea linealmente independiente y esas "fortalezas" estén en orden decreciente. Esta "fuerza" se captura exactamente en las diagonales de R correspondientes a este reordenamiento. De hecho, la fuerza es proporcional a este primer elemento. Lo que debe hacer es verificar qué diagonales de R en la versión reordenada son mayores que este primer coeficiente escalado por la tolerancia y las usa para determinar cuáles de las columnas correspondientes son linealmente independientes.

Sin embargo, voy a cambiar esto y determinar el punto en las diagonales R donde se encuentran las últimas columnas independientes posibles. Cualquier cosa después de este punto se consideraría linealmente dependiente . Esto es esencialmente lo mismo que verificar si las diagonales son menores que el umbral, pero estamos utilizando el reordenamiento de la matriz para nuestro beneficio.

En cualquier caso, poniendo lo que he mencionado en el código, esto es lo que debes hacer, suponiendo que tu matriz esté almacenada en A :

%// Step #1 - Define tolerance
tol = 1e-10;
%// Step #2 - Do QR Factorization
[Q, R, E] = qr(A,0);
diag_R = abs(diag(R)); %// Extract diagonals of R
%// Step #3 -
%// Find the LAST column in the re-arranged result that
%// satisfies the linearly independent property
r = find(diag_R >= tol*diag_R(1), 1, ''last'');
%// Step #4
%// Anything after r means that the columns are
%// linearly dependent, so let''s output those columns to the
%// user
idx = sort(E(r+1:end));

Tenga en cuenta que E será un vector de permutación, y supongo que quiere que estos sean ordenados, por eso los clasificamos después del punto en que los vectores ya no son linealmente independientes. Probemos esta teoría. Supongamos que tengo esta matriz:

A =
1 1 2 0
2 2 4 9
3 3 6 7
4 4 8 3

Puede ver que las dos primeras columnas son iguales, y la tercera columna es un múltiplo del primero o segundo. Simplemente tendrías que multiplicar uno por dos para obtener el resultado. Si ejecutamos el código anterior, esto es lo que obtengo:

idx =
1 2

Si también echas un vistazo a E , esto es lo que obtengo:

E =
4 3 2 1

Esto significa que la columna 4 fue la "mejor" columna independiente linealmente, seguida de la columna 3. Debido a que volvimos [1,2] como las columnas linealmente dependientes, esto significa que las columnas 1 y 2 tienen ambas [1,2,3,4] ya que sus columnas son un múltiplo escalar de alguna otra columna. En este caso, esta sería la columna 3 ya que las columnas 1 y 2 son la mitad de la columna 3.

¡Espero que esto ayude!

Método alternativo

Si no quiere hacer ninguna factorización QR , entonces le sugiero que reduzca su matriz a su forma Echelon reducida , y puede determinar los vectores de base que componen el espacio de columna de su matriz A Esencialmente, el espacio de columna le proporciona el conjunto mínimo de columnas que pueden generar todas las posibles combinaciones lineales de vectores de salida si tuviera que aplicar esta matriz mediante la multiplicación de matriz-vector. Puede determinar qué columnas forman el espacio de columna utilizando el comando rref . Proporcionaría una segunda salida a rref que le proporciona un vector de elementos que le indican qué columnas son linealmente independientes o forman una base del espacio de columnas para esa matriz. Como tal:

[B,RB] = rref(A);

RB le daría las ubicaciones de qué columnas para el espacio de columna y B sería la forma de escalón reducida de la matriz A Como desea encontrar aquellas columnas que son linealmente dependientes , le conviene devolver un conjunto de elementos que no contienen estas ubicaciones. Como tal, defina un vector linealmente creciente de 1 a tantas columnas como tenga, luego use RB para eliminar estas entradas en este vector y el resultado serían aquellas columnas linealmente dependientes que está buscando. En otras palabras:

[B,RB] = rref(A);
idx = 1 : size(A,2);
idx(RB) = [];

Al usar el código anterior, esto es lo que obtenemos:

idx =
2 3

Tenga en cuenta que identificamos que las columnas 2 y 3 son linealmente dependientes, lo que tiene sentido ya que ambas son múltiplos de la columna 1. La identificación de qué columnas son linealmente dependientes es diferente en comparación con el método de factorización QR, ya que QR ordena las columnas sobre la probabilidad de que esa columna en particular sea linealmente independiente. Debido a que las columnas 1 a 3 están relacionadas entre sí, no debería importar qué columna devuelve. Una de estas formas es la base de las otras dos.

No he probado la eficacia del uso de rref en comparación con el método QR. Sospecho que el rref elimina las filas gaussianas, donde la complejidad es peor en comparación con la factorización QR, ya que el algoritmo está altamente optimizado y es eficiente. Debido a que su matriz es bastante grande, me quedaré con el método QR, pero use rref todos modos y vea lo que obtiene.

Answer 2

Si normalizas cada columna dividiendo por su máximo, la proporcionalidad se convierte en igualdad. Esto hace que el problema sea más fácil.

Ahora, para probar la igualdad, puede usar un solo bucle (externo) sobre las columnas; el lazo interno se puede vectorizar fácilmente con bsxfun . Para una mayor velocidad, compare cada columna solo con las columnas a su derecha.

También para ahorrar algo de tiempo, la matriz de resultados se asigna previamente a un tamaño aproximado (debe establecer eso). Si el tamaño aproximado es incorrecto, la única penalización será un poco más lenta, pero el código funciona.

Como de costumbre, las pruebas de igualdad entre los valores de coma flotante deben incluir una tolerancia.

El resultado se da como una matriz de 2 columnas ( S ), donde cada fila contiene los índices de dos filas que son proporcionales.

A = [1 5 2 6 3 1
2 5 4 7 6 1
3 5 6 8 9 1]; %// example data matrix
tol = 1e-6; %// relative tolerance
A = bsxfun(@rdivide, A, max(A,[],1)); %// normalize A
C = size(A,2);
S = NaN(round(C^1.5),2); %// preallocate result to *approximate* size
used = 0; %// number of rows of S already used
for c = 1:C
ind = c+find(all(abs(bsxfun(@rdivide, A(:,c), A(:,c+1:end))-1)<tol));
u = numel(ind); %// number of columns proportional to column c
S(used+1:used+u,1) = c; %// fill in result
S(used+1:used+u,2) = ind; %// fill in result
used = used + u; %// update number of results
end
S = S(1:used,:); %// remove unused rows of S

En este ejemplo, el resultado es

S =
1 3
1 5
2 6
3 5

significa que la columna 1 es proporcional a la columna 3; la columna 1 es proporcional a la columna 5, etc.