Resolver $Ax=b$ en $L_1$ $|Ax-b|$ minimización

Question

Me gustaría resolver un sistema lineal $Ax = b$ bajo el $L_1$ restricción de la norma $\min(|Ax-b|)$ . Todo lo que puedo encontrar sobre $L_1$ La minimización es una forma de minimizar $|x|_1$ con sujeción a $Ax=b$ .

Quería utilizar la programación lineal en matlab para resolver este problema. Esto me llevó a resolver un nuevo sistema con linprog en matlab:

$$My = k$$

Así que hice algunas transformaciones, sabiendo que Linprog puede resolver el siguiente problema:

$$\min f(x) \ \ \textrm{s.t.} \ \ Ax \leq b$$

Quiero minimizar este problema:

$$\min ||My - k||_1$$

para $y$

Para convertir esto a la forma linprog, introduje un nuevo vector que quiero minimizar, $z$ :

$$\min \sum z$$

s.t.

$$My - k \leq z$$ $$-My + k \leq z$$

Sin embargo, las restricciones incluyen tanto la variable $y$ y las constantes $k$ . Podemos formularlo en la forma requerida creando un sistema mayor:

$Ax \leq b$ donde la matriz $A$ es:

$$\begin{pmatrix} M & -I \\ -M & I \end{pmatrix} $$

el vector $x$ es igual a $y$ y $z$ apilados: $$\begin{pmatrix} y \\ z \end{pmatrix} $$

y las restricciones son $k$ y $-k$ apilados: $$\begin{pmatrix} k \\ -k \end{pmatrix} $$

El problema es que esto no funciona. (Matlab no encuentra una buena solución). Me gustaría saber si hay otra forma de resolver este problema? ¿Si habéis oído hablar de algún código de Matlab que ya lo haga?

Muchas gracias

Answer 1

Si tienes min $\|Ax-b\|_1$ donde $A$ es $n\times m$ para $n>m$ y $b$ es $n\times 1$ (por lo que la suma de los valores absolutos de todos los componentes del vector resultado $y=Ax-b$ donde y es $n\times 1$ ) entonces puede resolver esto como un LP con linprog con un comando como este:

linprog([zeros(m,1);ones(n,1)],[+A,-eye(n);-A,-eye(n)],[b;-b])

que no es más que una traducción a Matlab del truco habitual descrito en Wikipedia .

He aquí un ejemplo de ajuste polinómico que compara las desviaciones mínimas absolutas con los mínimos cuadrados lineales.

clc
x = [0:1/100:1]';

%random polynomial of degree 4
poly = [7/960;-1/8;227/320;-23/15;1]/0.0583;
A = [x.^0,x.^1,x.^2,x.^3,x.^4];
b = A*poly;

% bigger k means less outliers. k = 2 has too many outliers
% for the least absolute deviations solution to lie on the initial
% polynomial
k = 3;
b(k*[1:101/k]) = x(k*[1:101/k]);

% linear least squares
poly = A\b;
y = A*poly;

% least absolute deviations
[n,m] = size(A);
poly = linprog([zeros(m,1);ones(n,1)],[+A,-eye(n);-A,-eye(n)],[b;-b]); poly = poly(1:m);
z = A*poly;

% blue lines: original data
% green lines: least squares solution
% magenta crosses: least absolute deviations solution
% red circles: outliers
plot(x,b,'b-',x,y,'g-',x,z,'m+',x(k*[1:101/k]),x(k*[1:101/k]),'ro')

% check least absolute deviations
check = fminsearch(@(x) norm(A*x-b,1), poly);
[[poly; norm(A*poly-b,1)], [check; norm(A*check-b,1)]]

Este es el aspecto de la salida

Answer 2

Te has olvidado de las restricciones de límite inferior $z \ge 0$ . Después de añadir esta restricción a LP linprog logró resolver su problema.