>> m=5;n=10;B=(1:5)';A=randn(m,n);x=zeros(n,1);x(B)=rand(m,1);b=A*x;c=randn(n,1); >> disp(b) -1.9935e-01 1.0275e+00 1.8235e+00 1.5965e+00 -2.0862e-01 >> disp([x c]) 4.0571e-01 -1.0181e+00 9.3547e-01 -1.8208e-01 9.1690e-01 1.5210e+00 4.1027e-01 -3.8439e-02 8.9365e-01 1.2274e+00 0 -6.9620e-01 0 7.5245e-03 0 -7.8289e-01 0 5.8694e-01 0 -2.5121e-01 >> disp(c'*x) 1.8924e+00 >> y=A(:,B)'\c(B); >> c-A'*y -6.6613e-16 4.4409e-16 2.2204e-16 4.5103e-16 0 -1.4265e-01 6.4881e-01 1.1772e+00 3.1167e+00 -1.1250e+00 >> s=n; >> T = A(:,B)\[A b] T = 1.0000e+00 0 5.9062e-17 0 0 3.3893e-01 1.0358e+00 1.2382e+00 3.7411e-02 -2.0118e+00 4.0571e-01 0 1.0000e+00 2.0179e-16 0 0 1.4004e+00 8.9248e-01 6.7111e+00 3.4227e+00 -6.2508e+00 9.3547e-01 0 0 1.0000e-00 0 0 -1.1701e-02 -1.2716e-01 -1.2765e+00 -1.3730e+00 5.3618e-01 9.1690e-01 0 0 2.4603e-17 1.0000e+00 0 -9.8239e-01 3.5974e-01 -1.3415e+00 -2.7095e+00 1.1087e+00 4.1027e-01 0 0 9.5787e-17 0 1.0000e+00 2.1624e-02 6.3791e-01 1.9655e+00 9.4228e-02 -2.5138e+00 8.9365e-01 >> [x(B) T(:,s) x(B)./T(:,s)] 4.0571e-01 -2.0118e+00 -2.0167e-01 9.3547e-01 -6.2508e+00 -1.4966e-01 9.1690e-01 5.3618e-01 1.7101e+00 4.1027e-01 1.1087e+00 3.7005e-01 8.9365e-01 -2.5138e+00 -3.5550e-01 >> r = 4; lammax=x(B(r))/T(r,s) lammax = 3.7005e-01 >> xnew1=zeros(n,1);xnew1(B)= x(B)-lammax*T(:,s);xnew1(s)=lammax xnew1 = 1.1502e+00 3.2485e+00 7.1849e-01 -5.5511e-17 1.8239e+00 0 0 0 0 3.7005e-01 >> Bnew1=[(1:3)';5;s] Bnew1 = 1 2 3 5 10 >> Tnew1=A(:,Bnew1)\A Tnew1 = 1.0000e+00 0 1.1627e-16 1.8145e+00 0 -1.4437e+00 1.6885e+00 -1.1960e+00 -4.8791e+00 -4.2150e-16 0 1.0000e+00 3.3794e-16 5.6379e+00 0 -4.1382e+00 2.9206e+00 -8.5230e-01 -1.1853e+01 -9.5321e-16 0 0 1.0000e-00 -4.8361e-01 0 4.6340e-01 -3.0113e-01 -6.2770e-01 -6.2609e-02 -1.1289e-16 0 0 1.5157e-16 2.2673e+00 1.0000e+00 -2.2058e+00 1.4536e+00 -1.0762e+00 -6.0491e+00 -4.2535e-16 0 0 2.2191e-17 9.0196e-01 0 -8.8607e-01 3.2447e-01 -1.2100e+00 -2.4438e+00 1.0000e-00 >> disp([c'*xnew1 c'*x]) 1.4761e+00 1.8924e+00 >> ynew1=A(:,Bnew1)'\c(Bnew1); >> c-A'*ynew1 2.2204e-16 -2.7756e-17 -2.2204e-16 1.0147e+00 0 -1.1395e+00 1.0138e+00 -1.8407e-01 3.6742e-01 1.6653e-16 >> s=6; >> disp([xnew1(Bnew1) Tnew1(:,s) xnew1(Bnew1)./Tnew1(:,s)]) 1.1502e+00 -1.4437e+00 -7.9670e-01 3.2485e+00 -4.1382e+00 -7.8501e-01 7.1849e-01 4.6340e-01 1.5505e+00 1.8239e+00 -2.2058e+00 -8.2686e-01 3.7005e-01 -8.8607e-01 -4.1763e-01 >> r=3; >> lammax1=xnew1(Bnew1(r))/Tnew1(r,s) lammax1 = 1.5505e+00 >> xnew2=zeros(n,1);xnew2(Bnew1) = xnew1(Bnew1)-lammax1*Tnew1(:,s);xnew2(s)=lammax1 xnew2 = 3.3885e+00 9.6648e+00 0 0 5.2439e+00 1.5505e+00 0 0 0 1.7439e+00 >> disp([c'*xnew2 c'*xnew1 c'*x]) -2.9068e-01 1.4761e+00 1.8924e+00 >> Bnew2=[Bnew1([1:r-1,r+1:5]);s] Bnew2 = 1 2 5 10 6 >> ynew2=A(:,Bnew2)'\c(Bnew2); >> c-A'*ynew2 ans = -2.2204e-16 5.5511e-17 2.4590e+00 -1.7449e-01 0 -1.1102e-16 2.7336e-01 -1.7276e+00 2.1346e-01 1.6653e-16 >> s=8; >> Tnew2=A(:,Bnew2)\A Tnew2 = 1.0000e+00 0 3.1154e+00 3.0790e-01 0 0 7.5038e-01 -3.1516e+00 -5.0741e+00 0 0 1.0000e+00 8.9302e+00 1.3192e+00 0 0 2.3150e-01 -6.4578e+00 -1.2412e+01 0 0 0 4.7600e+00 -3.4683e-02 1.0000e+00 0 2.0167e-02 -4.0641e+00 -6.3471e+00 0 0 0 1.9121e+00 -2.2775e-02 0 0 -2.5133e-01 -2.4102e+00 -2.5636e+00 1.0000e+00 0 0 2.1580e+00 -1.0436e+00 0 1.0000e+00 -6.4983e-01 -1.3546e+00 -1.3511e-01 0 >> disp([xnew2(Bnew2) Tnew2(:,s)]) 3.3885e+00 -3.1516e+00 9.6648e+00 -6.4578e+00 5.2439e+00 -4.0641e+00 1.7439e+00 -2.4102e+00 1.5505e+00 -1.3546e+00 OBJ UNBOUNDED FROM BELOW >> lam = 1e10; >> xnew3=zeros(n,1);xnew3(Bnew2)=xnew2(Bnew2)-lam*Tnew2(:,s);xnew3(s)=lam xnew3 = 3.1516e+10 6.4578e+10 0 0 4.0641e+10 1.3546e+10 0 1.0000e+10 0 2.4102e+10 >> disp([max(abs((A*xnew3-b))) norm(A*xnew3-b)/norm(A)/norm(xnew3)]); 2.5557e-05 8.2123e-17 >> disp([c'*xnew3 c'*xnew2 c'*xnew1 c'*x]) -1.7276e+10 -2.9068e-01 1.4761e+00 1.8924e+00