Universidad de ChileFacultad de Ciencias Fısicas y MatematicasDepartamento de Ingenierıa Matematica

Algoritmo del Pivote Complementario(o Metodo de Lemke)

Alumnos:Arnold GarcıaEmanuel Berrocal

Curso:MA-5701: Optimizacion No lineal

Profesores de Catedra:Jorge Amaya

Abderrahim Hantoute

Profesores Auxiliares:Francisco Santibanez

Carlos Flores G.

2 de julio de 2012

Indice

1. Introduccion y Motivacion 2

2. Problema Cuadratico 22.1. transformacion del problema cuadratico . . . . . . . . . . . . . . . . . . . . . . . . . . 3

3. Problema Complementario 43.1. Soluciones Complementarias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

4. Algoritmo del pivote complementario 5

5. Aplicaciones del Algoritmo 55.1. Ejemplo 1: Solucion Complementaria basica factible . . . . . . . . . . . . . . . . . . . 55.2. Ejemplo 2: Solucion Casi Complementaria basica factible . . . . . . . . . . . . . . . . 65.3. Propiedades del algoritmo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

6. Algunos resultados acerca del Metodo de Lemke. 76.1. Matrices copositivas y copositivas-plus. . . . . . . . . . . . . . . . . . . . . . . . . . . 76.2. Acerca de la convergencia del metodo de Lemke. . . . . . . . . . . . . . . . . . . . . . 8

7. Conclusiones 11

8. Anexo 128.1. Codigo matlab del algoritmo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128.2. Pruebas realizadas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

1. Introduccion y Motivacion

En este informe se estudiara una manera de resolver el problema cuadratico a traves de un pro-blema complementario a el, el cual surge de haber aplicado KKT a el problema cuadratico, dondeluego introduciremos el metodo de Lemke para resolver este problema y mostrar su convergencia enun numero finito de iteraciones, bajo supuestos adecuados, veremos que el algoritmo o bien para conuna solucion complementaria basica o bien concluiremos que el sistema original es inconsistente.

Es interesante notar como un problema que es cuadratico, se puede llevar en la practica a unproblema lineal, y mas aun, como aquellas restricciones de desigualdades que a veces son incomodasde trabajar, las podemos ver en el problema lineal (que llamaremos el problema complementario alproblema cuadratico), como restricciones de igualdad, luego de haber aplicado KKT.

Tambien lo interesante que veremos, es que el algoritmo que en una cantidad finita de pasos escapaz de resolver aquel problema complementario, lo hara usando una tecnica parecida a Simplexaplicado a un tableau, metodo que fue visto en el curso de Optimizacion.

2. Problema Cuadratico

A continuacion el problema cuadratico:

Definicion 2.1. Problema Cuadratico:

(P ) min1

2xtQx+ ctx

s.a Ax ≤ b

x ≥ 0.

Donde:Q: es una matriz simetrica de n x n.A: es una matriz de m x n.

Observacion 2.1. Si Q no es simetrica, entonces podemos definir Q = Q+Qt

2la cual es simetrica

y ademas xtQx = xtQx+xtQtx2

= 2xtQx2

= xtQx Por lo que siempre podremos considerar Q comosimetrica.

Dado que las restricciones son lineales podemos aplicar kkt a este problema.sea

L(−→x ,−→u ,−→v ) =1

2· xtQx+ ctx+ (Ax− b)u− xv

tenemos:

∇xL = ∇xL(−→x ,−→u ,−→v ) = Qx+ c+ Atu− v = 0

(ut, vt) ·(Ax− b−x

)= 0

2.1. transformacion del problema cuadratico

Como las restricciones son lineales por KKT, tendremos lo siguiente:

Qx+ c+ Atu− v = 0

ut(Ax− b) = 0

vtx = 0

u, v ≥ 0

⇔ Qx+ c+ Atu− v = 0

ut(Ax− b) = 0

vtx = 0

u, v ≥ 0

y = b− Axy ≥ 0

o bien

y + Ax = b

v − Atu−Qx = c

ytu = 0

vtx = 0

u, v ≥ 0

x, y ≥ 0

o bien (yv

)−[

0 −AAt Q

]·(ux

)=

(bc

)(∗)(

yv

)t

·(ux

)= 0(

yv

),

(ux

)≥ 0

3. Problema Complementario

A partir de lo anterior Definimos el siguiente problema:

Definicion 3.1.

(PC) w −Mz = q M : p × p (1)

wtz = 0 q : p × 1 (2)

w, z ≥ 0 w, z ∈ Rp (3)

Se dice que wj y zj para cada j = 1, . . . , p es un par de variables complementarias.

Observacion 3.1. para el problema anterior (∗), tendrıamos los siguiente:

w =

(yv

), q =

(bc

)z =

(ux

), M =

[0 −AAt Q

]Donde:x: Variable primaly: Holgura primalu: Variable dualv: Holgura dual

3.1. Soluciones Complementarias

Definicion 3.2. Una solucion de (1)-(2)-(3) se llama solucion complementaria basica factiblesi y solo si (w, z) es solucion basica factible de (1) a (3) y ademas una variable de cada par (wj, zj)es basica.

Observacion 3.2. Si q ≥ 0 entonces basta tomar w = q z = 0 y se tiene una solucion de (1), (2),(3). Con lo que ademas significa que x = 0, y = b, w = 0, v = c, es decir la solucion de (P ) es x=0.Si q � 0 entonces se agrega una nueva columna a la ecuacion (1):

w −Mz − z0e = q (4)

wjzj = 0 j = 1, . . . p (5)

wj, zj, z0 ≥ 0 j = 1, . . . p (6)

donde e = (1, . . . , 1)t

Ası este nuevo problema posee solucion a traves de iteraciones,pues tomamos z0 = max{−qj : 1 ≤ j ≤ z}, z = 0, y w = q + ez0, haciendo esto obtenemos unasolucion inicial para el sistema expuesto en (4), (5) y (6), las cuales son(solucion trivial):

z = 0

w = q + z0e

z0 = max{−qj : j = 1, . . . , p}La idea es partir de la solucion trivial y luego iterar hasta lograr una solucion de (4)-(5)-(6) conz0 = 0, ie, una solucion de (PC)

Definicion 3.3. Consideremos (4)-(5)-(6).Una solucion factible (w, z, z0) de este sistema se llama solucion basica factible casi complementariasi y solo si:

i) (w, z, z0) es solucion basica de (4), (6).

ii) Para algun s = 1, . . . , p, ws, zs son no basicas.

iii) z0 es basica

iv) Exactamente una variable del par (wj, zj) es basica, para todo j = 1, . . . , p, j 6= s.

4. Algoritmo del pivote complementario

Ya sabemos que si q ≥ 0 basta con tomar (w, z) = (q, 0) para que sea la solucion de (PC), porlo que esta algoritmo sirve para el caso cuando q � 0, por lo que tomamos el sistema de ecuacionesentre (4) y (5). Y se procede:

Definicion 4.1. Formamos el siguiente tableau[wI | −

z

M |−z01 | q

]Paso de Inicializacion: Sea qr = min{qi / i = 1, . . . , p} y pivotear en la columna z0 con fila r, parahacer entrar z0 a la base.Se define ys = zs y entramos al algoritmo:

Definicion 4.2.

1◦ paso Sea ds la columna actualizada en el tableau bajo la variable ys, si ds ≤ 0, ir al pa-so 4. En otro caso, determinar el ındice r para la minima proporcion siguiente: qr

drs=

min1≤i≤p

{qidis

: dis > 0}

donde q es la actualizacion de la columna del lado derecho. fi-

nalmente si la variable basica en la fila r es z0, ir al paso 3. En otro caso, ir al paso 2.

2◦ paso La variable basica en la fila r o es wl o zl, para algun l 6= s. La variable ys entra a la base yel tableau es actualizado pivoteando en la fila r y en la columna ys. Si la variable que solodeja la base es zl, entonces dejar ys = wl. Ir al paso 1

3◦ paso Aquı ys entra a la base, y z0 deja la base. Se pivotea en la columna ys y en la fila z0,produciendo una solucion complementaria basica factible. Entonces parar

4◦ paso Se encontro un rayo extremo. Un rayo R = {(w, z, z0) + λd / λ ≥ 0}, donde

d =

1, en la posicion de la variable que entraria(fila ys)−ds, en la fila de las actuales variables basicas0, en todos los demas

5. Aplicaciones del Algoritmo

A continuacion veremos los dos casos posible de solucion que nos puede dar el algoritmo del pivotecomplementario.

5.1. Ejemplo 1: Solucion Complementaria basica factible

Para este ejemplo, encontraremos una solucion complementaria basica factible. Tomaremos lossiguientes datos:

M =

0 1 11 2 11 1 1

q =

1−2

0

ahora empezamos la iteracion:

tendremos que en la 1◦ iteracion entrara z0 en la posicion r = 2, pues es el mınimo qr

variablesbasicas |w1 w2 w3| z1 z2 z3| z0| qw1 | 1 0 0| 0 −1 −1| −1| 1

w2 | 0 1 0|−1 −2 −1| −1 |−2 ← sale w2 y entra z0w3 | 0 0 1|−1 −1 −1| −1| 0

y por tanto saldra w2 de la base, quedando ys = z2, lo que ademas implica que ds = z2, luego eltableau, queda de la siguiente forma, ademas se indica lo que entrara y saldra de la base:

variablesbasicas |w1 w2 w3|z1 z2 z3|z0|qw1 | 1 −1 0| 1 1 0| 0|3z0 | 0 −1 0| 1 2 1| 1|2 ← sale z0 y entra z2w3 | 0 −1 1| 0 1 0| 0|2

Luego, nos queda el siguiente tableau:

variablesbasicas |w1 w2 w3|z1 z2 z3| z0|qw1 | 1 −1/2 0| 1 0 −1/2|−1/2|2z2 | 0 −1/2 0| 1 1 1/2| 1/2|1 ← Fin, pues z0 salio (estamosw3 | 0 −1/2 1| 0 0 −1/2|−1/2|1 en el paso 3◦ del algoritmo)

Con lo que se obtuvo la siguiente solucion complementaria basica factible:

(w, z)t = (2, 0, 1, 0, 1, 0)

Observacion 5.1. NOTAR: que despues de la 1◦ iteracion siempre tendremos que q ≥ 0

5.2. Ejemplo 2: Solucion Casi Complementaria basica factible

Para este ejemplo, encontraremos una solucion complementaria basica factible. Tomaremos lossiguientes datos:

M =

0 0 1 −10 0 −1 2−1 1 2 −21 −2 −2 2

q =

14−2−4

ahora empezamos la iteracion:

tendremos que en la 1◦ iteracion entrara z0, en la posicion r = 4, pues es el mınimo qr

variablesbasicas |w1 w2 w3 w4| z1 z2 z3 z4| z0| q|w1 | 1 0 0 0| 0 0 −1 1| −1| 1|w2 | 0 1 0 0| 0 0 1 −2| −1| 4|w3 | 0 0 1 0| 1 −1 −2 2| −1|−2|w4 | 0 0 0 1|−1 2 2 −2| −1 |−4| ← sale w4 y entra z0

y por tanto saldra w4 de la base, quedando ys = z4, lo que ademas implica que ds = z4, luego eltableau, queda de la siguiente forma, ademas se indica lo que entrara y saldra de la base:

variablesbasicas |w1 w2 w3 w4|z1 z2 z3 z4|z0|q|w1 | 1 0 0 −1| 1 −2 −3 3| 0|5|w2 | 0 1 0 −1| 1 −2 −1 0| 0|8|w3 | 0 0 1 −1| 2 −3 −4 4 | 0|2| ← sale w3 y entra z4z0 | 0 0 0 −1| 1 −2 −2 2| 1|4|

despues de esta iteracion tendremos que ys = z3, por lo que ds = z3, luego, el tableau, quedara:

variablesbasicas |w1 w2 w3 w4| z1 z2 z3 z4|z0| q|w1 | 1 0 −3/4 −1/4|−1/2 1/4 0 0| 0|7/2|← deberıa salir w1y deberıa entrar z3w2 | 0 1 0 −1| 1 −2 −1 0| 0| 8|z4 | 0 0 1/4 −1/4| 1/2 −3/4 −1 1| 0|1/2|z0 | 0 0 −1/2 −1/2| 0 −1/2 0 0| 1| 3|

Tal como sale en el tableau, deberıamos tener segun el algoritmo que z3 debiese salir, pero resulta,en este caso ds ≤ 0, por lo que estamos en el paso 4◦ del algoritmo, ası la solucion final es:

R = {(w, z, z0) = (7/2, 8, 0, 0, 0, 0, 0, 1/2, 3) + λ(0, 1, 0, 0, 0, 0, 0, 1, 1, 0) : λ ≥ 0}

Donde cada punto del conjunto R satisface (4)-(5)-(6).

5.3. Propiedades del algoritmo

Lema 5.1. Supongamos que cada solucion basica factible casi complementaria de (4)-(5)-(6) es nodegenerada, es decir, cada variable basica es positiva. Entonces los vectores (w, z, z0) generados porel algoritmo no se repiten, es decir el algoritmo se detiene al cabo de un numero finito de iteraciones.

Lema 5.2. Bajo las hipotesis de no degenerancia, supongamos que el algoritmo encuentra un rayoextremo, esto es:

(w, z, z0) + λ(w, z, z0)

con (w, z, z0) casi complementaria, basica factible.Entonces:

(i) (w, z, z0) 6= 0, w, z, z0 ≥ 0

(ii) w −Mz − z0e = 0

(iii) wtz = wtz = wtz = wtz = 0

(iv) z 6= 0

(v) zMz = −(etz)z0 ≤ 0

6. Algunos resultados acerca del Metodo de Lemke.

A continuacion y en la seccion siguiente se veran algunos resultados del metodo de Lemke, enparticular respecto a la convergencia de este metodo y la relacion entre las soluciones del algoritmoy las soluciones del problema complementario definido por las ecuaciones (1)− (2)− (3).Primero partiremos por caracterizar un tipo especial de matrices que usaremos en varios de los re-sultados que siguen a continuacion.

6.1. Matrices copositivas y copositivas-plus.

Partiremos con la siguiente definicion:

Definicion 6.1. Sea M matriz de p×p. Se dice que M es copositiva si y solo si ztMz > 0 ∀z > 0.Si ademas,

ztMz = 0

z > 0

}=⇒ (M t +M)z = 0

Se dice que M es copositiva-plus.

Recordemos que:

M =

[0 −AAt Q

]Ahora veamos en que casos la matriz M es copositiva y/o copositiva-plus.

Teorema 6.1. Sea Q simetrica, luego:

i) Q copositiva ⇒ M copositiva.

ii) Q copositiva-plus ⇒ M copositiva-plus.

Demostracion.

i) Sea zt = (xt, yt) > 0, luego

ztMz = (xt, yt)

[0 −AAt Q

](xy

)= (xt, yt)

(−Ay

Atx+Qy

)= ytQy > 0

ii) Supongamos que ztMz = 0 y z > 0, entonces

(M +M t)z =

[0 00 2Q

]z =

(0

2Qy

)Pero dado que

ztMz = 0 ⇒ ytQy = 0 y como Q es copositiva-plus

(Q+Qt)y = 0

⇒ 2Qy = 0

⇒ (M +M t)z = 0

Ademas se puede demostrar que:

i) Si Q es semidefinida positiva entonces M es copositiva-plus.

ii) Si Q = qij es tal que qij > 0 ∀i, j, entonces M es copositiva. Si ademas qii > 0 ∀i, entoncesM es copositiva-plus.

6.2. Acerca de la convergencia del metodo de Lemke.

Teorema 6.2. Consideremos el sistema (4)− (5)− (6) y supongamos que cada solucion complemen-taria basica factible es no degenerada. Supongamos ademas que M es copositiva-plus.Entonces el algoritmo de Lemke se detiene en una solucion complementaria basica factible del siste-ma (1)− (2)− (3) si y solo si el sistema (1)− (3) es consistente.

Demostracion. Supongamos que el algoritmo se detiene en una solucion complementaria basica fac-tible de (1)− (2)− (3). Entonces se encuentra (w, z, z0) pero como z0 esta fuera de la base se tieneque:

w −Mz = q

wtz = 0

w, z > 0

Por lo tanto se cumple que el sistema (1)− (3) es consistente.Ahora supongamos que el algoritmo se detiene en una solucion basica factible casi complementaria,luego la solucion de (4)− (5)− (6) es de la forma (w, z, z0) + λ(w, z, z0).Luego como estamos en un rayo extremo tenemos que z0 esta en la base y por la hipotesis de nodegenerancia se tiene que z0 > 0.

Ademas por el Lema (??)

z 6= 0, z > 0, ztMz = −(1tz)z0 6 0.

Y como M es copositiva entonces ztMz > 0 y luego (1tz)z0 = 0.Ası, dado que z 6= 0, entonces z0 = 0.Luego, como w −Mz − 1z0 = 0 tenemos w = Mz, lo que implica que Mz > 0.Ademas del hecho que ztMz = 0 y dado que M e copositiva-plus obtenemos que (M +M t)z = 0Es decir, M tz = −Mz 6 0, por lo tanto

M tz 6 0 (7)

Por otro lado, sabemos que:

w −Mz − 1z0 = q, luego

wz − ztM tz − (1tz)z0 = qtz

Del hecho que M tz = −Mz = −w se tiene:

wtz + ztw − (1tz)z0 = qtz

Y del Lema (??) se tiene wtz = ztw = 0, luego qtz = −(1tz)z0, ası

qtz < 0 (8)

De (7) y (8) y del hecho que z > 0 obtenemos

[−I,M t]z 6 0, −qtz > 0

Luego, del teorema de Farkas tenemos que:

[−I,M t]ty = −q, y > 0

No tiene solucion, es decir, notando y = (w, z)t el sistema

−w +Mz = −qw, z > 0

No tiene solucion, luego, (1)− (3) es inconsistente.

Teorema 6.3. Si M es copositiva-plus y todos sus elementos son positivos o nulos, salvo la diagonalen que todos son positivos, entonces el algoritmo alcanza una solucion complementaria basica factible.

Demostracion. tarea??????

Teorema 6.4. Consideremos el problema:

(PC) min1

2xtQx+ ctx = f(x)

s.a Ax ≤ b

x ≥ 0.

Y supongamos que existe al menos una solucion factible. Supongamos ademas que el metodo de Lemkees usado para resolver el problema complementario:

w −Mz = q

wtz = 0

w, z > 0

Si el problema complementario es no degenerado y se cumple al menos una de las 3 condicionessiguientes:

i) Q es semidefinida positiva y c = 0.

ii) Q es positiva definida.

iii) qij > 0 ∀i, j = 1, ..., n y q : ii > 0 ∀i = 1, ..., n

Entonces el algoritmo se detiene al cabo de un numero finito de iteraciones, en un punto de Kuhn-Tucker del problema (PC).Ademas, si Q es semidefinida positiva y el algoritmo se detiene en un rayo extremo de (PC), entonceses no acotado.

Demostracion.

i) Supongamos que el algoritmo termina en un rayo extremo. Como Q es semidefinida positiva en-tonces M es copositiva-plus, luego por el Teorema (6.1) se alcanza el rayo extremo solo si el sistemasiguiente no tiene solucion:

Ax+ y =b

−Qx− Atu+ v =0

x, y, u, v >0

⇔

{.w −Mz =q

w, z >0

Por el Teorema de Farkas, el siguiente sistema tiene solucion (d, e)

d 6 0 (9)

e 6 0 (10)

Ae > 0 (11)

Atd+Qe > 0 (12)

btd+ cte > 0 (13)

Multiplicando (12) por −et > 0 :

eAtd− etQe > 0

(Ae)td− etQe > 0⇒ etQe 6 0 pues (Ae)t 6 0

Como Q es semidefinida positiva etQe = 0, luego Q12 e = 0⇒ Qe = 0.

Como c = 0, de (13) btd > 0. Y dado que hemos supuesto que Ax 6 b, x > 0 tiene al menos unasolucion, sea x∗ > 0 esa solucion ,luego

Ax∗ + y∗ = b y∗ > 0

Tenemos

(Ax∗ + y∗)td > 0

x∗tAtd+ y∗td > 0

Pero de (12) se tiene que Atd 6 Qe = 0, luego

y∗td > 0 pues x∗ > 0

Lo que establece una contradiccion con (9).

ii) Razonando igual que en (i) obtenemos que e = 0

luego de (13) btd > 0⇔ x∗tAtd+ y∗td > 0

y de (12) Atd 6 0

Lo que implica la misma contradiccion.

iii) Razonando igual que en los caso anteriores resulta que etQe = 0 ⇒ e = 0 lo que da el mismoresultado que (ii).

7. Conclusiones

Hemos visto que ciertos problemas cuadraticos bajo, supuestos adecuados, es posible obtenersu problema complementario a traves de KKT, por ejemplo, alguna de esas ventajas fueron que:La matriz de restricciones fuera definida positiva o copositiva-plus, el algoritmo se detendra en unnumero finito de pasos.Ademas vimos las condiciones que se deben cumplir para que el problema complementario tengasolucion, haciendo al problema cuadratico un problema factible.

8. Anexo

En este se anexo se encontrara el codigo matlab utilizado para realizar el algoritmo del pivotecomplementario y ademas se encontraran algunas pruebas realizadas.

8.1. Codigo matlab del algoritmo

Aca se presenta el codigo Matlab para poder correr el algoritmo

%% Algoritmo de Pivote Complementario

function [w,z,z0,d]=pivote(M,q)

% M: Matriz de nxn.

% q: Vector en Rn.

n=length(M);

% base: si el elemento i de la base es k:

% k=1,...,n corresponde a los w_k

% k=n+1,...,2n corresponde a los z_k

% k=2n+1 corresponde a z0

base=1:n;

%%Paso 0 (Inicializacion)

if q>=0

w=q;

z=zeros(n,1);

z0=0;

d=zeros(2*n+1,1);

return

else

%Definir matrices

I=eye(n);

U=ones(n,1);

A=[I -M -U q];

[n m]=size(A);

%Elegir q_r

[minimo r]=min(q);

a=2*n+1;

%pivotear fila r, columna z0

A=Pivotear(A,r,a);

%variable que sale

var_s=base(r);

%Entra z0, actualizar base

base(r)=2*n+1;

end

while (true)

%%Paso 1

%Elegir variable que entra a la base

a=complementaria(var_s,n);

if A(:,a)<=0

%%% %%Paso 4

%salida

w=zeros(n,1);

z=zeros(n,1);

d=zeros(2*n+1,1);

for i=1:n

ind=base(i);

if ind<=n

w(ind)=A(i,m);

d(ind)=-A(i,a);

elseif ind==2*n+1

z0=A(i,m);

d(2*n+1)=-A(i,a);

else

z(ind-n)=A(i,m);

d(ind)=-A(i,a);

end

end

d(a)=1;

return

else

%Elegir fila que sale

r=radio_test(A,a);

%actualizar variable que sale

var_s=base(r);

if var_s==2*n+1

%%Paso 3

A=Pivotear(A,r,a);

base(r)=a;

%salida

z0=0;

d=zeros(2*n+1,1);

w=zeros(n,1);

z=zeros(n,1);

for i=1:n

ind=base(i);

if ind<=n

w(ind)=A(i,m);

else

z(ind-n)=A(i,m);

end

end

return

else

%%Paso 2

A=Pivotear(A,r,a);

base(r)=a;

end

end

end

end

%% Funcion para pivotear

function [B]=Pivotear(A,r,d)

[n m]=size(A);

P=eye(n);

P(:,r)=-A(:,d)/A(r,d);

P(r,r)=1/A(r,d);

B=P*A;

end

%% Funcion complementaria

function x=complementaria(y,n)

x=mod(y+n,2*n);

if x==0

x=2*n;

end

end

%% Funcion radio_test

function r=radio_test(A,a)

[n m]=size(A);

qq=A(:,m);

aa=A(:,a);

aa(aa<=0)=0;

[minimo r]=min(qq./aa);

end

8.2. Pruebas realizadas

%% Primer ejemplo

% Soluion complementaria basica factible:

% (w,z)’=(2 0 1 0 1 0)

clear

clc

M=[0 1 1;1 2 1;1 1 1];

q=[1;-2;0];

[w,z,z0,d]=pivote(M,q)

%% Segundo ejemplo

% Solucion basica factible casi complementaria

% (w,z,z0)’=(7/2 8 0 0 0 0 0 1/2 3) + L(0 1 0 0 0 0 1 1 0)

clear

clc

M=[0 0 1 -1;0 0 -1 2;-1 1 2 -2;1 -2 -2 2];

q=[1;4;-2;-4];

[w,z,z0,d]=pivote(M,q)