2012-11-24 07:24:58 +08:00
|
|
|
|
2012-12-12 03:14:36 +08:00
|
|
|
badscale = 1e-8;
|
2012-11-24 07:24:58 +08:00
|
|
|
|
|
|
|
|
% Create random Jacobian
|
2012-12-12 03:14:36 +08:00
|
|
|
% A = [
|
|
|
|
|
% badscale.*(rand(3)-.5) zeros(3) zeros(3,6)
|
|
|
|
|
% zeros(3) badscale.*(rand(3)-.5) zeros(3,6)
|
|
|
|
|
% rand(3,6)-.5 zeros(3,6)
|
|
|
|
|
% zeros(2,6) badscale.*(rand(2,6)-.5)
|
|
|
|
|
% eye(6) -eye(6)
|
|
|
|
|
% ];
|
|
|
|
|
% b = [
|
|
|
|
|
% repmat(badscale, 3,1)
|
|
|
|
|
% repmat(badscale, 3,1)
|
|
|
|
|
% repmat(1, 3,1)
|
|
|
|
|
% repmat(badscale, 2,1)
|
|
|
|
|
% repmat(1, 6,1)
|
|
|
|
|
% ];
|
|
|
|
|
|
|
|
|
|
Acoeffs = [
|
|
|
|
|
1 11 badscale
|
2016-04-12 03:11:29 +08:00
|
|
|
(1:10)' (1:10)' -Matrix::Ones(10,1)
|
|
|
|
|
(1:10)' (2:11)' Matrix::Ones(10,1)
|
2012-12-12 03:14:36 +08:00
|
|
|
]';
|
|
|
|
|
A = full(sparse(Acoeffs(1,:), Acoeffs(2,:), Acoeffs(3,:)));
|
|
|
|
|
b = zeros(size(A,1), 1);
|
|
|
|
|
|
|
|
|
|
% A = zeros(2*15, 2*10);
|
|
|
|
|
% for j = 2:10
|
|
|
|
|
% A(j*2-3 : j*2-2, j*2-3:j*2) = rand(2,4) - 0.5;
|
|
|
|
|
% end
|
|
|
|
|
% A(19:20, 1:2) = badscale * (rand(2) - 0.5);
|
|
|
|
|
% A(21:22, [ 3 4 11 12 ]) = rand(2,4) - 0.5;
|
|
|
|
|
% A(23:24, [ 1 2 5 6 ]) = rand(2,4) - 0.5;
|
|
|
|
|
% A(25:26, [ 5 6 17 18 ]) = rand(2,4) - 0.5;
|
|
|
|
|
% A(27:28, [ 9 10 13 14 ]) = rand(2,4) - 0.5;
|
|
|
|
|
% A(29:30, [ 7 8 15 16 ]) = rand(2,4) - 0.5;
|
|
|
|
|
% b = zeros(2*15, 1);
|
|
|
|
|
% b(1:18) = rand(18,1) - 0.5;
|
|
|
|
|
% b(19:20) = badscale * (rand(2,1) - 0.5);
|
|
|
|
|
% b(21:30) = rand(10,1) - 0.5;
|
|
|
|
|
|
|
|
|
|
% Acoeffs = [
|
|
|
|
|
% 1 1 1 1 1 1 2 2 2 2 2 3 3 3 3 4 4 4 5 5 6 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3 4 4 4 4 4 4 5 5 5 5 5 5 6 6 6 6 6 6 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3 4 4 4 4 4 4 5 5 5 5 5 5 6 6 6 6 6 6 1 2 3 4 5 6 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3 4 4 4 4 4 4 5 5 5 5 5 5 6 6 6 6 6 6 1 2 3 4 5 6 1 2 3 4 5 6 7 7 7 8 8 9 7 7 7 7 7 7 8 8 8 8 8 8 9 9 9 9 9 9 7 7 7 7 7 7 8 8 8 8 8 8 9 9 9 9 9 9 7 8 9 7 7 7 7 7 7 8 8 8 8 8 8 9 9 9 9 9 9 7 8 9 7 8 9 10 10 10 10 10 10 11 11 11 11 11 12 12 12 12 13 13 13 14 14 15 10 10 10 10 10 10 11 11 11 11 11 11 12 12 12 12 12 12 13 13 13 13 13 13 14 14 14 14 14 14 15 15 15 15 15 15 10 11 12 13 14 15 10 10 10 10 10 10 11 11 11 11 11 11 12 12 12 12 12 1
|
|
|
|
|
% 39 40 41 42 43 44 40 41 42 43 44 41 42 43 44 42 43 44 42 44 42 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 10 11 12 13 14 15 10 11 12 13 14 15 10 11 12 13 14 15 10 11 12 13 14 15 10 11 12 13 14 15 10 11 12 13 14 15 4 5 6 7 8 9 4 5 6 7 8 9 4 5 6 7 8 9 4 5 6 7 8 9 4 5 6 7 8 9 4 5 6 7 8 9 38 38 38 38 38 38 32 33 34 35 36 37 32 33 34 35 36 37 32 33 34 35 36 37 32 33 34 35 36 37 32 33 34 35 36 37 32 33 34 35 36 37 31 31 31 31 31 31 45 45 45 45 45 45 1 2 3 1 3 1 10 11 12 13 14 15 10 11 12 13 14 15 10 11 12 13 14 15 4 5 6 7 8 9 4 5 6 7 8 9 4 5 6 7 8 9 38 38 38 32 33 34 35 36 37 32 33 34 35 36 37 32 33 34 35 36 37 31 31 31 45 45 45 10 11 12 13 14 15 10 12 13 14 15 12 13 14 15 13 14 15 13 14 14 4 5 6 7 8 9 4 5 6 7 8 9 4 5 6 7 8 9 4 5 6 7 8 9 4 5 6 7 8 9 4 5 6 7 8 9 38 38 38 38 38 38 32 33 34 35 36 37 32 33 34 35 36 37 32 33 34 35 36 3
|
|
|
|
|
% 2287.81 -4.4988e-05 0.00396445 -0.00100866 6.88479 -0.382977 2287.81 -0.0282353 -6.88394 -0.000338518 -0.210359 2287.43 0.384701 0.211465 -0.000303552 0.409179 574.421 -0.33448 -0.086675 573.083 572.434 0.000175554 -1.36517 0.05849 1.36507 0.000112365 0.0189909 -0.058652 -0.0188071 3.14326e-05 -0.0247666 -71.2266 0.0198345 0.00567252 -0.0219283 -71.2443 -71.2327 0.02605 -0.00547148 42.1211 427.311 -13.0572 0.0738133 -0.0064156 0.0071859 425.771 -40.7137 40.8349 -0.00719115 -0.0660657 0.0058001 -39.2512 16.9899 430.388 -0.00195469 -0.00782391 3.41702e-05 0.532335 -0.69246 -1.01117 -0.51164 0.0871608 0.000338874 -0.0645723 -1.86503 0.0118343 -0.184545 0.00037516 -0.353041 -5.38585 0.256705 2.13342 -0.0527851 -0.910563 0.00685038 -2288.11 -0.319849 1.49778 0.00131762 -6.68744 0.436231 0.318124 -2288.11 1.17676 6.68538 -0.000179114 0.288793 -1.53826 -1.20334 -2289.06 -0.435568 -0.295666 0.000650458 -8.20257 -0.00495529 -1.40144 -1.19995 -721.834 0.999144 0.0115191 1.39146 0.00206421 0.241407 0.578931 -715.253 -0.00426102 8.23871 0.0319042 -712.673 -0.687971 0.134907 1.78921e-05 -0.000506192 -0.00661338 0.00136493 0.000168206 -0.00104316 0.819305 0.0028821 0.00347505 -0.000315376 -0.197366 -0.0521121 0.00209998 0.830532 0.0283978 0.198642 0.000553189 -0.0787171 0.00964912 0.0871257 1.43774 0.0514481 0.0844689 -0.000378355 -3.5015 -0.389049 -0.260307 0.791191 147.416 -0.664512 -5.2737 0.465659 -0.142763 -0.154666 -0.578578 142.184 0.51651 7.02367 3.46957 140.242 0.689716 -0.134755 0.00138855 -0.00029297 -4.60677e-07 -5.85457e-06 0.00341196 -1.84581e-07 -0.169162 -0.143293 -0.103757 0.971775 -9.55431 -17.9962 0.0254247 37.4517 -0.0209152 -0.00499967 37.3707 37.3282 1.62285 10.577 -1.26238 -10.5926 1.04948 -0.975285 -0.167226 -2.22364 0.016771 0.755791 -0.415351 -10.4962 -16.5666 1.33474 1.44369 -1.04654 -11.0613 0.329368 -101.978 -0.0119957 -0.19074 -0.671635 -143.392 0.564983 1.84747 0.281797 0.00213787 0.13017 0.483793 -138.977 0.0625116 102.371 -1.91922 -137.339 -0.576069 0.111915 0.00133071 0.000165007 -0.00101693 -3.36239 -0.383114 -0.255928 0.672042 143.394 -0.564846 -5.15865 0.455748 -0.139196 -0.130105 -0.483446 138.993 0.507372 6.83154 3.4007 137.342 0.577777 -0.111759 3.88677e-05 0.00503247 9.4483e-06 2.00827 0.481391 -16.6988 1.36255 2075.39 -21.4475 18.6028 0.118129 0.928177 2065.47 16.8895 -0.0624046 -21.3263 0.313145 1618.3 -0.696466 1.66511 -0.11786 -2.34839 1.47139 121.564 83.1783 -0.597982 73.8588 453.099 -43.4199 16.1018 -0.292859 -1.68581 0.851333 44.5456 454.444 -37.9494 4.16356 -0.404957 0.0926619 -10.0948 45.2476 595.572 -1.03824 0.501659 0.00710507 -3.69312 -0.816079 0.57251 0.184137 -0.868655 -11.6822 -105.54 10.8294 1.39689 -1.49873 -15.7957 2.6502 12.1566 136.446 -25.4427 -11.7228 1.02809 -0.247507 0.0527668 -0.0277195 2.42197 0.00384665 0.0175 -0.11796 -230.563 20.5743 -15.8678 0.290267 1.68717 -0.781925 -21.9179 -231.021 30.1778 -4.17005 0.402651 -0.090854 1.80197 -56.0322 -523.372 0.824811 -0.599641 0.0020112
|
|
|
|
|
% ];
|
|
|
|
|
% Ab = sparse(Acoeffs(1,:), Acoeffs(2,:), Acoeffs(3,:));
|
|
|
|
|
% A = full(Ab(:,1:end-1));
|
|
|
|
|
% b = full(Ab(:,end));
|
2012-11-24 07:24:58 +08:00
|
|
|
|
|
|
|
|
% Solve with MATLAB
|
|
|
|
|
x_true = A \ b;
|
|
|
|
|
|
|
|
|
|
% Form Hessian
|
|
|
|
|
L = A'*A;
|
|
|
|
|
l = A'*b;
|
|
|
|
|
|
|
|
|
|
% Eliminate badly-conditioned system. Use this custom choleskyNaive function
|
|
|
|
|
% that does no internal pivoting or scaling to simulate the case where the
|
|
|
|
|
% bad conditioning occurs between fronts.
|
|
|
|
|
R = choleskyNaive(L);
|
|
|
|
|
if isempty(R)
|
|
|
|
|
fprintf('Cholesky failed on badly-scaled system with a negative pivot\n');
|
|
|
|
|
else
|
|
|
|
|
clear opts
|
|
|
|
|
opts.LT = true;
|
|
|
|
|
d = linsolve(R', l, opts);
|
|
|
|
|
|
|
|
|
|
% Solve badly-conditioned system.
|
|
|
|
|
clear opts
|
|
|
|
|
opts.UT = true;
|
|
|
|
|
x = linsolve(R, d, opts);
|
|
|
|
|
|
|
|
|
|
% Compute error
|
|
|
|
|
fprintf('Solution error with Cholesky on badly-scaled system: %g\n', ...
|
|
|
|
|
max(abs(x - x_true)));
|
|
|
|
|
end
|
|
|
|
|
|
|
|
|
|
% Form scaled system from Jacobian
|
2012-12-12 03:14:36 +08:00
|
|
|
% Aqr = [
|
|
|
|
|
% A(1:9, :)
|
|
|
|
|
% -eye(6) eye(6) ];
|
|
|
|
|
% Si = full(qr(sparse(Aqr)));
|
|
|
|
|
% Si = Si(1:size(Si,2), :);
|
|
|
|
|
% S = eye(size(Si,1)) / Si;
|
|
|
|
|
|
|
|
|
|
S = eye(size(A,2));
|
|
|
|
|
|
|
|
|
|
% Aqr = [
|
|
|
|
|
% A(19:20, :)
|
|
|
|
|
% -eye(6) eye(6) ];
|
|
|
|
|
|
|
|
|
|
% Aqr = A(1:20, :);
|
|
|
|
|
% Si = full(qr(sparse(Aqr)));
|
|
|
|
|
% Si = Si(1:size(Si,2), :);
|
|
|
|
|
% S = eye(size(Si,1)) / Si;
|
|
|
|
|
|
|
|
|
|
% bandwidth = 6;
|
|
|
|
|
% Lchol = diag(diag(L));
|
|
|
|
|
% for i = 1:bandwidth
|
|
|
|
|
% Lchol = Lchol + diag(diag(L,i),i) + diag(diag(L,-i),-i);
|
|
|
|
|
% end
|
|
|
|
|
% Si = chol(Lchol);
|
|
|
|
|
% S = inv(Si);
|
|
|
|
|
|
2012-11-24 07:24:58 +08:00
|
|
|
%Si = diag(sqrt(diag(L)));
|
2012-12-12 03:14:36 +08:00
|
|
|
% S = [
|
|
|
|
|
% eye(6) zeros(6)
|
|
|
|
|
% eye(6) eye(6) ];
|
2012-11-24 07:24:58 +08:00
|
|
|
Ap = A * S;
|
|
|
|
|
Lp = Ap'*Ap;
|
|
|
|
|
lp = Ap'*b;
|
|
|
|
|
|
|
|
|
|
% Eliminate scaled system
|
|
|
|
|
Rp = choleskyNaive(Lp);
|
|
|
|
|
if isempty(Rp)
|
|
|
|
|
fprintf('Cholesky failed on scaled system with a negative pivot\n');
|
|
|
|
|
else
|
|
|
|
|
clear opts
|
|
|
|
|
opts.LT = true;
|
|
|
|
|
dp = linsolve(Rp', lp, opts);
|
|
|
|
|
|
|
|
|
|
% Solve rescaled system.
|
|
|
|
|
clear opts
|
|
|
|
|
opts.UT = true;
|
|
|
|
|
xp = S * linsolve(Rp, dp, opts);
|
|
|
|
|
|
|
|
|
|
% Compute error
|
|
|
|
|
fprintf('Solution error with Cholesky from rescaled Jacobian: %g\n', ...
|
|
|
|
|
max(abs(xp - x_true)));
|
|
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
% Form scaled system from Hessian
|
|
|
|
|
Lph = S' * (L * S);
|
|
|
|
|
lph = S' * l;
|
|
|
|
|
|
|
|
|
|
% Eliminate scaled system
|
|
|
|
|
Rph = choleskyNaive(Lph);
|
|
|
|
|
if isempty(Rph)
|
|
|
|
|
fprintf('Cholesky failed on scaled system with a negative pivot\n');
|
|
|
|
|
else
|
|
|
|
|
clear opts
|
|
|
|
|
opts.LT = true;
|
|
|
|
|
dph = linsolve(Rph', lph, opts);
|
|
|
|
|
|
|
|
|
|
% Solve rescaled system.
|
|
|
|
|
clear opts
|
|
|
|
|
opts.UT = true;
|
|
|
|
|
xph = S * linsolve(Rph, dph, opts);
|
|
|
|
|
|
|
|
|
|
% Compute error
|
|
|
|
|
fprintf('Solution error with Cholesky from rescaled Hessian: %g\n', ...
|
|
|
|
|
max(abs(xph - x_true)));
|
|
|
|
|
end
|
|
|
|
|
|
