performance - subplot matlab
MATLAB Optimización de la ortogonalización ponderada de Gram-Schmidt (3)
Tengo una función en MATLAB que realiza la ortogonalización Gram-Schmidt con una ponderación muy importante aplicada a los productos internos (no creo que la función incorporada de MATLAB lo admita). Esta función funciona bien por lo que puedo decir, sin embargo, es demasiado lenta en matrices grandes. ¿Cuál sería la mejor manera de mejorar esto?
He intentado convertir a un archivo MEX pero pierdo la paralelización con el compilador que estoy usando y, por lo tanto, es más lento.
Estaba pensando en ejecutarlo en una GPU ya que las multiplicaciones de elementos están altamente paralelizadas. (Pero prefiero que la implementación sea fácilmente transportable)
¿Puede alguien vectorizar este código o hacerlo más rápido? No estoy seguro de cómo hacerlo elegantemente ...
Sé que las mentes de stackoverflow aquí son asombrosas, considera esto como un desafío :)
Función
function [Q, R] = Gram_Schmidt(A, w)
[m, n] = size(A);
Q = complex(zeros(m, n));
R = complex(zeros(n, n));
v = zeros(n, 1);
for j = 1:n
v = A(:,j);
for i = 1:j-1
R(i,j) = sum( v .* conj( Q(:,i) ) .* w ) / ...
sum( Q(:,i) .* conj( Q(:,i) ) .* w );
v = v - R(i,j) * Q(:,i);
end
R(j,j) = norm(v);
Q(:,j) = v / R(j,j);
end
end
donde A es una matriz mxn de números complejos y w es un vector mx 1 de números reales.
Embotellamiento
Esta es la expresión para R(i,j) que es la parte más lenta de la función (no es 100% seguro si la notación es correcta):
donde w es una función de peso no negativo. El producto interno ponderado se menciona en varias páginas de Wikipedia, esta es una sobre la función de peso y esta sobre las funciones ortogonales .
Reproducción
Puedes producir resultados usando el siguiente script:
A = complex( rand(360000,100), rand(360000,100));
w = rand(360000, 1);
[Q, R] = Gram_Schmidt(A, w);
donde A y w son las entradas.
Velocidad y computación
Si utiliza la secuencia de comandos anterior, obtendrá los resultados del generador de perfiles a lo siguiente:
Resultado de la prueba
Puede probar los resultados comparando una función con la anterior utilizando el siguiente script:
A = complex( rand( 100, 10), rand( 100, 10));
w = rand( 100, 1);
[Q , R ] = Gram_Schmidt( A, w);
[Q2, R2] = Gram_Schmidt2( A, w);
zeros1 = norm( Q - Q2 );
zeros2 = norm( R - R2 );
donde Gram_Schmidt es la función descrita anteriormente y Gram_Schmidt2 es una función alternativa. Los resultados zeros1 y zeros2 deberían estar muy cerca de cero.
Nota:
Intenté acelerar el cálculo de R(i,j) con lo siguiente pero en vano ...
R(i,j) = ( w'' * ( v .* conj( Q(:,i) ) ) ) / ...
( w'' * ( Q(:,i) .* conj( Q(:,i) ) ) );
1)
Mi primer intento de vectorización:
function [Q, R] = Gram_Schmidt1(A, w)
[m, n] = size(A);
Q = complex(zeros(m, n));
R = complex(zeros(n, n));
for j = 1:n
v = A(:,j);
QQ = Q(:,1:j-1);
QQ = bsxfun(@rdivide, bsxfun(@times, w, conj(QQ)), w.'' * abs(QQ).^2);
for i = 1:j-1
R(i,j) = (v.'' * QQ(:,i));
v = v - R(i,j) * Q(:,i);
end
R(j,j) = norm(v);
Q(:,j) = v / R(j,j);
end
end
Desafortunadamente, resultó ser más lento que la función original.
2)
Entonces me di cuenta de que las columnas de esta matriz intermedia QQ se construyen de forma incremental y que las anteriores no se modifican. Así que aquí está mi segundo intento:
function [Q, R] = Gram_Schmidt2(A, w)
[m, n] = size(A);
Q = complex(zeros(m, n));
R = complex(zeros(n, n));
QQ = complex(zeros(m, n-1));
for j = 1:n
if j>1
qj = Q(:,j-1);
QQ(:,j-1) = (conj(qj) .* w) ./ (w.'' * (qj.*conj(qj)));
end
v = A(:,j);
for i = 1:j-1
R(i,j) = (v.'' * QQ(:,i));
v = v - R(i,j) * Q(:,i);
end
R(j,j) = norm(v);
Q(:,j) = v / R(j,j);
end
end
Técnicamente no se realizó ninguna vectorización importante; Solo he precalculado los resultados intermedios y moví el cálculo fuera del bucle interno.
Basado en un rápido punto de referencia, esta nueva versión es definitivamente más rápida:
% some random data
>> M = 10000; N = 100;
>> A = complex(rand(M,N), rand(M,N));
>> w = rand(M,1);
% time
>> timeit(@() Gram_Schmidt(A,w), 2) % original version
ans =
1.2444
>> timeit(@() Gram_Schmidt1(A,w), 2) % first attempt (vectorized)
ans =
2.0990
>> timeit(@() Gram_Schmidt2(A,w), 2) % final version
ans =
0.4698
% check results
>> [Q,R] = Gram_Schmidt(A,w);
>> [Q2,R2] = Gram_Schmidt2(A,w);
>> norm(Q-Q2)
ans =
4.2796e-14
>> norm(R-R2)
ans =
1.7782e-12
EDITAR:
Luego de los comentarios, podemos reescribir la segunda solución para deshacernos de if-statmenet, al mover esa parte al final del bucle externo (es decir, inmediatamente después de calcular la nueva columna Q(:,j) , calculamos y almacenamos la QQ(:,j) correspondiente QQ(:,j) ).
La función es idéntica en la salida, y el tiempo tampoco es tan diferente; ¡El código es un poco más corto!
function [Q, R] = Gram_Schmidt3(A, w)
[m, n] = size(A);
Q = zeros(m, n, ''like'',A);
R = zeros(n, n, ''like'',A);
QQ = zeros(m, n, ''like'',A);
for j = 1:n
v = A(:,j);
for i = 1:j-1
R(i,j) = (v.'' * QQ(:,i));
v = v - R(i,j) * Q(:,i);
end
R(j,j) = norm(v);
Q(:,j) = v / R(j,j);
QQ(:,j) = (conj(Q(:,j)) .* w) ./ (w.'' * (Q(:,j).*conj(Q(:,j))));
end
end
Tenga en cuenta que utilicé la sintaxis de zeros(..., ''like'',A) (nuevo en las últimas versiones de MATLAB). Esto nos permite ejecutar la función sin modificar en la GPU (suponiendo que tenga la Caja de herramientas de computación paralela):
% CPU
[Q3,R3] = Gram_Schmidt3(A, w);
contra
% GPU
AA = gpuArray(A);
[Q3,R3] = Gram_Schmidt3(AA, w);
Lamentablemente en mi caso, no fue más rápido. De hecho, fue muchas veces más lento correr en la GPU que en la CPU, pero valió la pena un tiro :)
Es posible vectorizar esto por lo que solo se necesita un bucle. El cambio fundamental importante del algoritmo original es que si intercambias los bucles interno y externo, puedes vectorizar la proyección del vector de referencia a todos los vectores restantes. Trabajando con la solución de @Amro , encontré que un bucle interno es en realidad más rápido que la resta de matrices. No entiendo por qué esto sería. Sincronizando esto con la solución de @Amro , es aproximadamente un 45% más rápido.
function [Q, R] = Gram_Schmidt5(A, w)
Q = A;
n_dimensions = size(A, 2);
R = zeros(n_dimensions);
R(1, 1) = norm(Q(:, 1));
Q(:, 1) = Q(:, 1) ./ R(1, 1);
for i = 2 : n_dimensions
Qw = (Q(:, i - 1) .* w)'' * Q(:, (i - 1) : end);
R(i - 1, i : end) = Qw(2:end) / Qw(1);
%% Surprisingly this loop beats the matrix multiply
for j = i : n_dimensions
Q(:, j) = Q(:, j) - Q(:, i - 1) * R(i - 1, j);
end
%% This multiply is slower than above
% Q(:, i : end) = ...
% Q(:, i : end) - ...
% Q(:, i - 1) * R(i - 1, i : end);
R(i, i) = norm(Q(:,i));
Q(:, i) = Q(:, i) ./ R(i, i);
end
Hay una larga discusión aquí, pero, para saltar a la respuesta. Ha ponderado el numerador y el denominador del cálculo de R mediante un vector w. La ponderación se produce en el bucle interno y consiste en un producto de punto triple, un punto Q punto w en el numerador y Q punto Q punto w en el denominador. Si realiza un cambio, creo que el código se ejecutará significativamente más rápido. Escriba num = (A punto sqrt (w)) punto (Q punto sqrt (w)) y escriba den = (Q punto sqrt (w)) punto (Q punto sqrt (w)). Eso mueve los cálculos del producto (A dot sqrt (w)) y (Q dot sqrt (w)) fuera del bucle interno.
Me gustaría dar una descripción de la formulación de la ortogonalización de Gram Schmidt, que, con suerte, además de brindar una solución computacional alternativa, brinda una perspectiva adicional de la ventaja de la OSG.
Los "objetivos" de la OSG son dos. Primero, para habilitar la solución de una ecuación como Ax = y, donde A tiene muchas más filas que columnas. Esta situación ocurre con frecuencia cuando se miden datos, ya que es fácil medir más datos que el número de estados. El enfoque para el primer objetivo es reescribir A como QR, de modo que las columnas de Q sean ortogonales y normalizadas, y R sea una matriz triangular. El algoritmo que proporcionaste, creo, logra el primer objetivo. La Q representa el espacio base de la matriz A, y R representa la amplitud de cada espacio base requerido para generar cada columna de A.
El segundo objetivo de la OSG es clasificar los vectores de base en orden de importancia. Este es el paso que no has hecho. Y, mientras se incluye este paso, puede aumentar el tiempo de solución, los resultados identificarán qué elementos de x son importantes, de acuerdo con los datos contenidos en las mediciones representadas por A.
Pero creo que con esta implementación, la solución es más rápida que el enfoque que presentó.
Aij = Qij Rij donde Qj son ortonormales y Rij es triangular superior, Ri, j> i = 0. Qj son los vectores de base ortogonal para A, y Rij es la participación de cada Qj para crear una columna en A. Entonces,
A_j1 = Q_j1 * R_1,1
A_j2 = Q_j1 * R_1,2 + Q_j2 * R_2,2
A_j3 = Q_j1 * R_1,3 + Q_j2 * R_2,3 + Q_j3 * R_3,3
Por inspección, puede escribir
A_j1 = ( A_j1 / | A_j1 | ) * | A_j1 | = Q_j1 * R_1,1
Luego proyectas Q_j1 en cada columna A para obtener los elementos R_1, j
R_1,2 = Q_j1 dot Aj2
R_1,3 = Q_j1 dot Aj3
...
R_1,j(j>1) = A_j dot Q_j1
Luego restas los elementos del proyecto de Q_j1 de las columnas de A (esto establecería la primera columna en cero, por lo que puedes ignorar la primera columna).
for j = 2,n
A_j = A_j - R_1,j * Q_j1
end
Ahora se ha eliminado una columna de A, se determinó el primer vector de base ortonormal, Q, j1, y se determinó la contribución del primer vector de base a cada columna, R_1, j, y la contribución del primer vector de base tiene Se ha restado de cada columna. Repita este proceso para las columnas restantes de A para obtener las columnas restantes de Q y las filas de R.
for i = 1,n
R_ii = |A_i| A_i is the ith column of A, |A_i| is magnitude of A_i
Q_i = A_i / R_ii Q_i is the ith column of Q
for j = i, n
R_ij = | A_j dot Q_i |
A_j = A_j - R_ij * Q_i
end
end
Estás tratando de ponderar las filas de A, con w. Aquí hay un enfoque. Yo normalizaría w, e incorporaría el efecto en R. Usted "eliminó" los efectos de w multiplicando y dividiendo por w. Una alternativa a "eliminar" el efecto es normalizar la amplitud de w a uno.
w = w / | w |
for i = 1,n
R_ii = |A_i inner product w| # A_i inner product w = A_i .* w
Q_i = A_i / R_ii
for j = i, n
R_ij = | (A_i inner product w) dot Q_i | # A dot B = A'' * B
A_j = A_j - R_ij * Q_i
end
end
Otro enfoque para implementar w es normalizar w y luego premultuar cada columna de A por w. Eso pesa limpiamente las filas de A y reduce el número de multiplicaciones.
Usar lo siguiente puede ayudar a acelerar tu código
A inner product B = A .* B
A dot w = A'' w
(A B)'' = B''A''
A'' conj(A) = |A|^2
Lo anterior se puede vectorizar fácilmente en matlab, más o menos como está escrito.
Pero, le falta la segunda parte de la clasificación de A, que le dice qué estados (elementos de x en A x = y) están representados significativamente en los datos
El procedimiento de clasificación es fácil de describir, pero le dejaré trabajar los detalles de la programación. El procedimiento anterior esencialmente asume que las columnas de A están en orden de importancia, y la primera columna se resta de todas las columnas restantes, luego la 2ª columna se resta de las columnas restantes, etc. La primera fila de R representa la contribución de primera columna de Q a cada columna de A. Si sumas el valor absoluto de la primera fila de contribuciones de R, obtienes una medida de la contribución de la primera columna de Q a la matriz A. Entonces, solo evalúas cada columna de Q. A como la primera (o siguiente) columna de Q, y determine el puntaje de clasificación de la contribución de esa columna de Q a las columnas restantes de A. Luego, seleccione la columna A que tenga el rango más alto como la siguiente columna de Q. La codificación de esto se reduce esencialmente a la estimación previa de la siguiente fila de R, para cada columna restante en A, para determinar qué magnitud R clasificada tiene la mayor amplitud. Tener un vector índice que represente el orden de la columna original de A será beneficioso. Al clasificar los vectores de base, terminas con los vectores de base "principales" que representan A, que suele ser mucho más pequeño que el número de columnas en A.
Además, si clasifica las columnas, no es necesario calcular todas las columnas de R. Cuando sepa qué columnas de A no contienen información útil, no hay un beneficio real para mantener esas columnas.
En dinámica estructural, un enfoque para reducir el número de grados de libertad es calcular los valores propios, asumiendo que usted tiene valores representativos para la matriz de masa y rigidez. Si lo piensa, el enfoque anterior se puede usar para "calcular" las matrices M y K (y C) a partir de la respuesta medida, y también para identificar las "formas de respuesta de medición" que están representadas significativamente en los datos. Estas son diferentes, y potencialmente más importantes, que las formas de modo. Por lo tanto, puede resolver problemas muy difíciles, es decir, la estimación de matrices de estado y el número de grados de libertad representados, a partir de la respuesta medida, mediante el enfoque anterior. Si leyó sobre N4SID, él hizo algo similar, excepto que usó SVD en lugar de GSO. No me gusta la descripción técnica de N4SID, demasiado énfasis en la notación de proyección vectorial, que es simplemente un producto puntual.
Puede haber uno o dos errores en la información anterior, estoy escribiendo esto desde lo alto de mi cabeza, antes de irme a trabajar. Entonces, verifique el algoritmo / ecuaciones a medida que implementa ... Buena suerte
Volviendo a su pregunta, de cómo optimizar el algoritmo cuando pesa con w. Aquí hay un algoritmo GSO básico, sin la clasificación, escrito compatible con su función.
Tenga en cuenta que el siguiente código está en octava, no en matlab. Hay algunas diferencias menores.
function [Q, R] = Gram_Schmidt_2(A, w)
[m, n] = size(A);
Q = complex(zeros(m, n));
R = complex(zeros(n, n));
# Outer loop identifies the basis vectors
for j = 1:n
aCol = A(:,j);
# Subtract off the basis vector
for i = 1:(j-1)
R(i,j) = ctranspose(Q(:,j)) * aCol;
aCol = aCol - R(i,j) * Q(:,j);
end
amp_A_col = norm(aCol);
R(j,j) = amp_A_col;
Q(:,j) = aCol / amp_A_col;
end
end
Para obtener su algoritmo, solo cambie una línea. Pero, pierdes mucha velocidad porque "ctranspose (Q (:, j)) * aCol" es una operación vectorial pero "sum (aCol. * Conj (Q (:, i)). * W)" es una fila operación.
function [Q, R] = Gram_Schmidt_2(A, w)
[m, n] = size(A);
Q = complex(zeros(m, n));
R = complex(zeros(n, n));
# Outer loop identifies the basis vectors
for j = 1:n
aCol = A(:,j);
# Subtract off the basis vector
for i = 1:(j-1)
# R(i,j) = ctranspose(Q(:,j)) * aCol;
R(i,j) = sum( aCol .* conj( Q(:,i) ) .* w ) / ...
sum( Q(:,i) .* conj( Q(:,i) ) .* w );
aCol = aCol - R(i,j) * Q(:,j);
end
amp_A_col = norm(aCol);
R(j,j) = amp_A_col;
Q(:,j) = aCol / amp_A_col;
end
end
Puede volver a cambiarlo a una operación vectorial al ponderar aCol y Q mediante el cuadrado de w.
function [Q, R] = Gram_Schmidt_3(A, w)
[m, n] = size(A);
Q = complex(zeros(m, n));
R = complex(zeros(n, n));
Q_sw = complex(zeros(m, n));
sw = w .^ 0.5;
for j = 1:n
aCol = A(:,j);
aCol_sw = aCol .* sw;
# Subtract off the basis vector
for i = 1:(j-1)
# R(i,j) = ctranspose(Q(:,i)) * aCol;
numTerm = ctranspose( Q_sw(:,i) ) * aCol_sw;
denTerm = ctranspose( Q_sw(:,i) ) * Q_sw(:,i);
R(i,j) = numTerm / denTerm;
aCol_sw = aCol_sw - R(i,j) * Q_sw(:,i);
end
aCol = aCol_sw ./ sw;
amp_A_col = norm(aCol);
R(j,j) = amp_A_col;
Q(:,j) = aCol / amp_A_col;
Q_sw(:,j) = Q(:,j) .* sw;
end
end
Como lo señaló JacobD, la función anterior no se ejecuta más rápido. Posiblemente tome tiempo para crear matrices adicionales. Otra estrategia de agrupación para el producto triple es agrupar w con conj (Q). Espero que esto sea más rápido ...
function [Q, R] = Gram_Schmidt_4(A, w)
[m, n] = size(A);
Q = complex(zeros(m, n));
R = complex(zeros(n, n));
for j = 1:n
aCol = A(:,j);
for i = 1:(j-1)
cqw = conj(Q(:,i)) .* w;
R(i,j) = ( transpose( aCol ) * cqw) ...
/ (transpose( Q(:,i) ) * cqw);
aCol = aCol - R(i,j) * Q(:,i);
end
amp_A_col = norm(aCol);
R(j,j) = amp_A_col;
Q(:,j) = aCol / amp_A_col;
end
end
Aquí hay una función de controlador para sincronizar diferentes versiones.
function Gram_Schmidt_tester_2
nSamples = 360000;
nMeas = 100;
nMeas = 15;
A = complex( rand(nSamples,nMeas), rand(nSamples,nMeas));
w = rand(nSamples, 1);
profile on;
[Q1, R1] = Gram_Schmidt_basic(A);
profile off;
data1 = profile ("info");
tData1=data1.FunctionTable(1).TotalTime;
approx_zero1 = A - Q1 * R1;
max_value1 = max(max(abs(approx_zero1)));
profile on;
[Q2, R2] = Gram_Schmidt_w_Orig(A, w);
profile off;
data2 = profile ("info");
tData2=data2.FunctionTable(1).TotalTime;
approx_zero2 = A - Q2 * R2;
max_value2 = max(max(abs(approx_zero2)));
sw=w.^0.5;
profile on;
[Q3, R3] = Gram_Schmidt_sqrt_w(A, w);
profile off;
data3 = profile ("info");
tData3=data3.FunctionTable(1).TotalTime;
approx_zero3 = A - Q3 * R3;
max_value3 = max(max(abs(approx_zero3)));
profile on;
[Q4, R4] = Gram_Schmidt_4(A, w);
profile off;
data4 = profile ("info");
tData4=data4.FunctionTable(1).TotalTime;
approx_zero4 = A - Q4 * R4;
max_value4 = max(max(abs(approx_zero4)));
profile on;
[Q5, R5] = Gram_Schmidt_5(A, w);
profile off;
data5 = profile ("info");
tData5=data5.FunctionTable(1).TotalTime;
approx_zero5 = A - Q5 * R5;
max_value5 = max(max(abs(approx_zero5)));
profile on;
[Q2a, R2a] = Gram_Schmidt2a(A, w);
profile off;
data2a = profile ("info");
tData2a=data2a.FunctionTable(1).TotalTime;
approx_zero2a = A - Q2a * R2a;
max_value2a = max(max(abs(approx_zero2a)));
profshow (data1, 6);
profshow (data2, 6);
profshow (data3, 6);
profshow (data4, 6);
profshow (data5, 6);
profshow (data2a, 6);
sprintf(''Time for %s is %5.3f sec with %d samples and %d meas, max value is %g'',
data1.FunctionTable(1).FunctionName,
data1.FunctionTable(1).TotalTime,
nSamples, nMeas, max_value1)
sprintf(''Time for %s is %5.3f sec with %d samples and %d meas, max value is %g'',
data2.FunctionTable(1).FunctionName,
data2.FunctionTable(1).TotalTime,
nSamples, nMeas, max_value2)
sprintf(''Time for %s is %5.3f sec with %d samples and %d meas, max value is %g'',
data3.FunctionTable(1).FunctionName,
data3.FunctionTable(1).TotalTime,
nSamples, nMeas, max_value3)
sprintf(''Time for %s is %5.3f sec with %d samples and %d meas, max value is %g'',
data4.FunctionTable(1).FunctionName,
data4.FunctionTable(1).TotalTime,
nSamples, nMeas, max_value4)
sprintf(''Time for %s is %5.3f sec with %d samples and %d meas, max value is %g'',
data5.FunctionTable(1).FunctionName,
data5.FunctionTable(1).TotalTime,
nSamples, nMeas, max_value5)
sprintf(''Time for %s is %5.3f sec with %d samples and %d meas, max value is %g'',
data2a.FunctionTable(1).FunctionName,
data2a.FunctionTable(1).TotalTime,
nSamples, nMeas, max_value2a)
end
En mi vieja computadora portátil, en Octave, los resultados son
ans = Time for Gram_Schmidt_basic is 0.889 sec with 360000 samples and 15 meas, max value is 1.57009e-16
ans = Time for Gram_Schmidt_w_Orig is 0.952 sec with 360000 samples and 15 meas, max value is 6.36717e-16
ans = Time for Gram_Schmidt_sqrt_w is 0.390 sec with 360000 samples and 15 meas, max value is 6.47366e-16
ans = Time for Gram_Schmidt_4 is 0.452 sec with 360000 samples and 15 meas, max value is 6.47366e-16
ans = Time for Gram_Schmidt_5 is 2.636 sec with 360000 samples and 15 meas, max value is 6.47366e-16
ans = Time for Gram_Schmidt2a is 0.905 sec with 360000 samples and 15 meas, max value is 6.68443e-16
Estos resultados indican que el algoritmo más rápido es el algoritmo sqrt_w arriba a 0,39 segundos, seguido de la agrupación de conj (Q) con w (arriba) a 0,452 segundos, luego la versión 2 de Amro Solution a 0,905 segundos, luego el algoritmo original en la pregunta en 0.952, luego una versión 5 que intercambia filas / columnas para ver si se presenta el almacenamiento de la fila (código no incluido) en 2.636 seg. Estos resultados indican que la división sqrt (w) entre A y Q es la solución más rápida. Pero estos resultados no son consistentes con el comentario de JacobD acerca de que sqrt (w) no es más rápido.