Problema inverso de las pdes.

Question

Lineal inversa problema está dada por:

$\ \mathbf{d}=\mathbf{A}\mathbf{m}+\mathbf{e}$

donde d: datos observados, Una: la teoría de operador, m: desconocido modelo y e: error.

Para minimizar el efecto del ruido; un Error Mínimo Cuadrado (LSE) modelo de estimación se utiliza comúnmente:

$\ \mathbf{\tilde{m}}=(\mathbf{A^\top A})^{-1}\mathbf{A^\top d}$

Como un ejemplo de problema que estoy teniendo en cuenta el modelo:

$\ T(x,z) = 0.5*sin(2\pi*x) - z$

Mis datos observados es el ruido de las derivadas parciales de T:

$\ P(x,z)=T_x+e$

$\ Q(x,z)=T_z+e$

Mediante la introducción de la matriz de diferencias finitas operadores puedo formular estas dos ecuaciones en derivadas parciales lineal de ecuaciones:

$\ \mathbf{P=D_x*T} $

$\ \mathbf{Q=D_z*T} $

Abajo está la LSE solución de estas ecuaciones:

Hay un gran problema con esta solución: en el modelo de la línea de donde T=0 pasa entre z=-0.5 y z=0.5, mientras que en la estimación va entre z=-1 y z=1. Algunos "ampliación" de la coordenada z parece estar teniendo lugar.

¿Por qué está sucediendo esto?

Aparte de esto la solución es bastante bueno, salvo por el bloque patrón. Este patrón parece estar relacionado con el muestreo; el más fino es el muestreo de la más fina sea la forma de bloques patrón.

¿Cuál es la explicación para este patrón?

Con el fin de obtener una suave solución que he probado con un suavizado de regularización de la restricción: $\ a*\mathbf{L*T}=0$

donde a es un parámetro de escala y L es la diferencia finita de Laplace de la matriz.

Me di cuenta de un trade-off entre la suavidad y la corrección según el parámetro de escala. Para los grandes parámetros de la escala de la onda sinusoidal patrón de la modelo se convierte en una onda cuadrada patrón:

¿Cómo puedo determinar la idoneidad de una?

Los problemas inversos son a menudo perturbada por las inestabilidades que hacen necesaria la creación de un amortiguador de regularización plazo. Este no parece ser el caso para este problema. Sin embargo, ya que Una no es cuadrada no podía encontrar a su condición de número con cond (a) en Matlab.

Hay alguna manera en Matlab para encontrar si una LSE problema con un rectangulares de la matriz es inestable?

Por CIERTO, Aquí está el código de Matlab he utilizado:

% number of samples in x and z direction:
n = 101;

% amount of noise added to partial derivatives:
e = 0.5;

x = linspace(-1,1,n);
z = linspace(-1,1,n);
dx = x(2) - x(1);
dz = z(2) - z(1);
[X,Z] = meshgrid(x,z);

% Hidden model:
subplot(1,3,1);
T = 0.5*sin(2*pi*X) - Z;
imagesc(x,z,T);
xlabel('x');
ylabel('z');
title('T = 0.5*sin(2*pi*x) - z');
colorbar;

% Observed noisy data:
P = pi*cos(2*pi*X) + e*randn(n,n); % Tx
Q = -ones(n,n)     + e*randn(n,n); % Tz

subplot(1,3,2);
imagesc(x,z,P);
xlabel('x');
ylabel('z');
title('P = Tx + e');
colorbar;
subplot(1,3,3);
imagesc(x,z,Q);
xlabel('x');
ylabel('z');
title('Q = Tz + e');
colorbar;

% Central difference matrices:
E = sparse(2:n, 1:n-1, 1, n, n);
D = (E' - E) / 2;
I = speye(n); 
Dx = kron(D,I)/dx;
Dz = kron(I,D)/dz;

A = []; b = [];
A = [A; Dx]; b = [b; P(:)]; 
A = [A; Dz]; b = [b; Q(:)];

% Least Squares Solution:
figure;
subplot(1,2,1);
imagesc(x,z,T);
xlabel('x');
ylabel('z');
title('T');
colorbar;

tlse = lsqr(A, b, 1e-3, 1000*100);
Tlse = reshape(tlse, [n n]);
subplot(1,2,2);
imagesc(x,z,Tlse);
xlabel('x');
ylabel('z');
title('Tlse');
colorbar;

figure;
% Least Squares Solution regularized by smoothing Laplace operator:
D2 = E + E' + sparse(1:n, 1:n, -2, n, n);
Dxx = kron(D2,I)/(dx^2);
Dzz = kron(I,D2)/(dz^2);
L = Dxx + Dzz;
ns = 3;
si = 1;
for si = 1 : ns;
    subplot(1,ns,si);
    % regularization:
    a = (si - 1)*5e-4;
    Ar = [A; a*L]; br = [b; zeros(n^2,1)];
    tlse = lsqr(Ar, br, 1e-3, 1000*100);
    Tlse = reshape(tlse, [n n]);
    imagesc(x,z,Tlse);
    str = sprintf('Tlse a=%g',a);
    title(str);
    si = si + 1;
end

EDITAR: Mi central a diferencia de la matriz D tenían un borde problema. Como un ejemplo de un tamaño de n=5: fue:

D =

         0    0.5000         0         0         0
   -0.5000         0    0.5000         0         0
         0   -0.5000         0    0.5000         0
         0         0   -0.5000         0    0.5000
         0         0         0   -0.5000         0

Reescribió este a

D =

   -1.0000    1.0000         0         0         0
   -0.5000         0    0.5000         0         0
         0   -0.5000         0    0.5000         0
         0         0   -0.5000         0    0.5000
         0         0         0   -1.0000    1.0000

Con más correcta, la diferencia de las matrices puedo obtener mucho mejores resultados:

Y esto es sin ningún tipo de regularización. Lo siento por el descuido. He aprendido mucho más difícil :-). Gracias por sus respuestas!

Answer 1

Factor de 2

No puedo estar seguro, pero parece que sus centrales diferencia de la matriz tiene un factor adicional de 1/2 que no debería estar allí. Que podrían provocar que la estimación incorrecta de la escala por (1/2*1/2)^-1*(1/2)=2.

La condición

La condición es la condición de número de $A^T A$, ya que es la matriz que se está invertida.

La elección del parámetro de regularización

La elección del parámetro de regularización es un problema clásico, para el que existe una vasta literatura y muchos métodos. Un buen libro sobre el tema es "la Regularización de los Problemas Inversos" por Engl, Hanke, y Neubauer (aunque carece de la evolución reciente). Uno particularmente simple pero eficaz método es la discrepancia principio.

Discrepancia principio: Elija el parámetro de regularización de modo que el error mínimo cuadrado residual es de aproximadamente el mismo tamaño que el nivel de ruido. Por ejemplo, si el nivel de ruido es $\delta$, comenzar con $k=0$ y resolver el problema inverso con $\alpha=1/2^k$. A continuación, compruebe si $||Am-d||<\delta$. Si es así, parar, si no lo es, establecer $k=k+1$ y repetir.

¿Cuál es el origen de la forma de bloques de ruido

Para mal planteado problemas inversos, las funciones propias del operador son generalmente más altos y mayor frecuencia oscilatoria funciones, con más y más pequeños autovalores decae a cero. Cuando intenta invertir, los de alta frecuencia de los modos de obtener amplificado.

Si usted tiene un vistazo a su ruido, verá parece que alguien ha añadido una alta frecuencia de la sinusoide en la parte superior de la verdadera función - hubo un poco de ruido en alta frecuencia de la componente senoidal, y fue amplificado enormemente por $(A^TA)^{-1}$. También hubo ruido en los otros componentes, pero no se amplifica como mucho, ya que sus autovalores son mayores (modo 1/autovalor es menor).

Si desea una demostración visual de lo que está pasando, calcular eig($A^TA$), y la trama de la última eigenfunction.

Malla de la regularización

Una gruesa malla puede actuar como la regularización, si la naturaleza de la malla se promedian los componentes de frecuencia alta. Mallas en realidad puede ser elegido para el propósito de regularización - hacerlo de forma adaptativa en realidad es un área de investigación actual!

Answer 2

También me gustaría considerar si sus aproximaciones a las primeras derivadas no introducen un desacoplamiento impar a través de la diferenciación central que utiliza: cambie la diferenciación a una aproximación unilateral para ver si eso mejora su solución no regularizada.