% gdama09
% W. Menke (2024), Geophysical Data Analysis and Inverse Theory
% with MATLAB(R) or Python, 5th Edition
% Code History:
% 2022-10-01, Imported from 4th Edition by W. Menke
% 2023-04-06, Edit and integration with text
% 2023-07-11, Added Householder example
% gdama09_01
% depicts a 3D box containing a vector
clear all;
fprintf('gdama09_01\n');
gdama09_01
% independent variable x
Nx = 51;
xmin = 0;
xmax = 1;
Dx = (xmax-xmin)/(Nx-1);
x = xmin + Dx*[0:Nx-1]';
% independent variable x
Ny = 51;
ymin = 0;
ymax = 1;
Dy = (ymax-ymin)/(Ny-1);
y = ymin + Dy*[0:Ny-1]';
% independent variable x
Nz = 51;
zmin = 0;
zmax = 1;
Dz = (zmax-zmin)/(Nz-1);
z = zmin + Dz*[0:Nz-1]';
% make grid
[XX, YY, ZZ] = meshgrid( x, y, z );
% model parameter graph
% make a "ball" in 3-space by contouring
% a Normal distribution centered on the
% center of the ball
% parameters for Normal distribution
rbar = [0.7, 0.8, 0.9]';
sd=0.1;
C = (sd^2)*eye(3,3);
CI = inv(C);
DC = det(C);
norm = ( ((2*pi)^(3/2)) * sqrt(DC) );
% normal distribution
PP = zeros(Nx,Ny,Nz);
for i=(1:Nx)
for j=(1:Ny)
for k=(1:Nz)
r = [XX(i,j,k), YY(i,j,k), ZZ(i,j,k)]';
PP(i,j,k) = exp(-0.5*(r-rbar)'*CI*(r-rbar))/norm;
end
end
end
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [xmin, xmax, ymin, ymax, zmin, zmax]');
% improvise outline of 3D box
plot3( [xmin,xmin], [ymin,ymin], [zmin,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymin,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmax], [ymin,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymin,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymin,ymin], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
xlabel('m_1');
ylabel('m_2');
zlabel('m_3');
% plot a line representing a vector
maxP=max(max(max(PP)));
plot3( [0, rbar(1)], [0, rbar(2)], [0, rbar(3)], 'k-', 'LineWidth', 3 );
% pretty crazy way to draw an arrowhead!
tangent = rbar/sqrt(rbar'*rbar);
per1 = cross( tangent, [0, 0, 1]' );
per1 = per1/sqrt(per1'*per1);
per2 = cross( tangent, per1 );
per2 = per2/sqrt(per2'*per2);
L = 0.05;
v1 = rbar - L*tangent + 0.25*L*per1;
v2 = rbar - L*tangent - 0.25*L*per1;
v3 = rbar - L*tangent + 0.25*L*per2;
v4 = rbar - L*tangent - 0.25*L*per2;
plot3( [rbar(1), v1(1)], [rbar(2), v1(2)], [rbar(3), v1(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v2(1)], [rbar(2), v2(2)], [rbar(3), v2(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v3(1)], [rbar(2), v3(2)], [rbar(3), v3(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v4(1)], [rbar(2), v4(2)], [rbar(3), v4(3)], 'k-', 'LineWidth', 3 );
% contour the ball
p=patch(isosurface( XX, YY, ZZ, PP, 0.95*maxP ));
isonormals(XX,YY,ZZ,PP, p)
set(p, 'FaceColor', 'red', 'FaceAlpha', 0.5, 'EdgeColor', 'none');
% set view, lighting, etc
daspect([1 1 1])
view(3)
camlight; lighting phong
fprintf('Caption: Model parameters in the space of possible models, depicted as a\n');
Caption: Model parameters in the space of possible models, depicted as a
fprintf('point (red) and vector (black)\n');
point (red) and vector (black)
% data graph
% parameters for Normal distribution
rbar = [0.9, 0.6, 0.7]';
sd=0.1;
C = (sd^2)*eye(3,3);
CI = inv(C);
DC = det(C);
norm = ( ((2*pi)^(3/2)) * sqrt(DC) );
% Normal distribution
PP = zeros(Nx,Ny,Nz);
for i=(1:Nx)
for j=(1:Ny)
for k=(1:Nz)
r = [XX(i,j,k), YY(i,j,k), ZZ(i,j,k)]';
PP(i,j,k) = exp(-0.5*(r-rbar)'*CI*(r-rbar))/norm;
end
end
end
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [xmin, xmax, ymin, ymax, zmin, zmax]');
% improvise outline of 3D box
plot3( [xmin,xmin], [ymin,ymin], [zmin,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymin,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmax], [ymin,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymin,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymin,ymin], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
xlabel('d_1');
ylabel('d_2');
zlabel('d_3');
% draw a vector
maxP=max(max(max(PP)));
plot3( [0, rbar(1)], [0, rbar(2)], [0, rbar(3)], 'k-', 'LineWidth', 3 );
% pretty crazy way to draw an arrowhead!
tangent = rbar/sqrt(rbar'*rbar);
per1 = cross( tangent, [0, 0, 1]' );
per1 = per1/sqrt(per1'*per1);
per2 = cross( tangent, per1 );
per2 = per2/sqrt(per2'*per2);
L = 0.05;
v1 = rbar - L*tangent + 0.25*L*per1;
v2 = rbar - L*tangent - 0.25*L*per1;
v3 = rbar - L*tangent + 0.25*L*per2;
v4 = rbar - L*tangent - 0.25*L*per2;
plot3( [rbar(1), v1(1)], [rbar(2), v1(2)], [rbar(3), v1(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v2(1)], [rbar(2), v2(2)], [rbar(3), v2(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v3(1)], [rbar(2), v3(2)], [rbar(3), v3(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v4(1)], [rbar(2), v4(2)], [rbar(3), v4(3)], 'k-', 'LineWidth', 3 );
% contour the ball
p=patch(isosurface( XX, YY, ZZ, PP, 0.95*maxP ));
isonormals(XX,YY,ZZ,PP, p)
set(p, 'FaceColor', 'blue', 'FaceAlpha', 0.5, 'EdgeColor', 'none');
% set view angle, lighting, etc.
daspect([1 1 1])
view(3)
camlight; lighting phong
fprintf('Caption: Datum in the space of possible measurements, depicted as a\n');
Caption: Datum in the space of possible measurements, depicted as a
fprintf('point (blue) and vector (black)\n');
point (blue) and vector (black)
% gdama09_02
% draw vectors in 3D space
clear all;
fprintf('gdama09_02\n');
gdama09_02
% independent variable, x
xmin = 0;
xmax = 1;
% independent variable, y
ymin = 0;
ymax = 1;
% independent variable, z
zmin = 0;
zmax = 1;
% graph 1
% these three vectors span 3D space
m1 = [0.7, 0.8, 0.8]';
m2 = [0.5, 0.7, 0.8]';
m3 = [0.8, 0.6, 0.7]';
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [xmin, xmax, ymin, ymax, zmin, zmax]');
% improvise outline of 3D box
plot3( [xmin,xmin], [ymin,ymin], [zmin,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymin,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmax], [ymin,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymin,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymin,ymin], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
xlabel('m_1');
ylabel('m_2');
zlabel('m_3');
% pretty crazy way to draw arrowheads!
% first arrowhead
rbar=m1; % point of the arrow
plot3( [0, rbar(1)], [0, rbar(2)], [0, rbar(3)], 'k-', 'LineWidth', 3 );
tangent = rbar/sqrt(rbar'*rbar); % tangent to the vector
per1 = cross( tangent, [0, 0, 1]' ); % a vector perpendicular to tangent
per1 = per1/sqrt(per1'*per1);
per2 = cross( tangent, per1 ); % another vector perpendicular to tangent
per2 = per2/sqrt(per2'*per2);
% arrowhead has 4 lines, each sloping back from the tip
L = 0.05;
v1 = rbar - L*tangent + 0.25*L*per1;
v2 = rbar - L*tangent - 0.25*L*per1;
v3 = rbar - L*tangent + 0.25*L*per2;
v4 = rbar - L*tangent - 0.25*L*per2;
plot3( [rbar(1), v1(1)], [rbar(2), v1(2)], [rbar(3), v1(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v2(1)], [rbar(2), v2(2)], [rbar(3), v2(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v3(1)], [rbar(2), v3(2)], [rbar(3), v3(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v4(1)], [rbar(2), v4(2)], [rbar(3), v4(3)], 'k-', 'LineWidth', 3 );
% second arrowhead
rbar=m2;
plot3( [0, rbar(1)], [0, rbar(2)], [0, rbar(3)], 'k-', 'LineWidth', 3 );
tangent = rbar/sqrt(rbar'*rbar);
per1 = cross( tangent, [0, 0, 1]' );
per1 = per1/sqrt(per1'*per1);
per2 = cross( tangent, per1 );
per2 = per2/sqrt(per2'*per2);
L = 0.05;
v1 = rbar - L*tangent + 0.25*L*per1;
v2 = rbar - L*tangent - 0.25*L*per1;
v3 = rbar - L*tangent + 0.25*L*per2;
v4 = rbar - L*tangent - 0.25*L*per2;
plot3( [rbar(1), v1(1)], [rbar(2), v1(2)], [rbar(3), v1(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v2(1)], [rbar(2), v2(2)], [rbar(3), v2(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v3(1)], [rbar(2), v3(2)], [rbar(3), v3(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v4(1)], [rbar(2), v4(2)], [rbar(3), v4(3)], 'k-', 'LineWidth', 3 );
% third arrowhead
rbar=m3;
plot3( [0, rbar(1)], [0, rbar(2)], [0, rbar(3)], 'k-', 'LineWidth', 3 );
tangent = rbar/sqrt(rbar'*rbar);
per1 = cross( tangent, [0, 0, 1]' );
per1 = per1/sqrt(per1'*per1);
per2 = cross( tangent, per1 );
per2 = per2/sqrt(per2'*per2);
L = 0.05;
v1 = rbar - L*tangent + 0.25*L*per1;
v2 = rbar - L*tangent - 0.25*L*per1;
v3 = rbar - L*tangent + 0.25*L*per2;
v4 = rbar - L*tangent - 0.25*L*per2;
plot3( [rbar(1), v1(1)], [rbar(2), v1(2)], [rbar(3), v1(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v2(1)], [rbar(2), v2(2)], [rbar(3), v2(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v3(1)], [rbar(2), v3(2)], [rbar(3), v3(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v4(1)], [rbar(2), v4(2)], [rbar(3), v4(3)], 'k-', 'LineWidth', 3 );
% draw triangle connecting vector tips
plot3( [m1(1), m2(1)] , [m1(2), m2(2)], [m1(3), m2(3)], 'r:', 'LineWidth', 2 );
plot3( [m2(1), m3(1)] , [m2(2), m3(2)], [m2(3), m3(3)], 'r:', 'LineWidth', 2 );
plot3( [m3(1), m1(1)] , [m3(2), m1(2)], [m3(3), m1(3)], 'r:', 'LineWidth', 2 );
% set view
view(3);
fprintf('Caption: The three vectors (black) span the 3D space\n');
Caption: The three vectors (black) span the 3D space
% graph 2
% these three vectors dont span 3D space
m1 = [0.7, 0.8, 0.9]';
m2 = [0.5, 0.7, 0.8]';
m3 = 0.5*(m1+m2);
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [xmin, xmax, ymin, ymax, zmin, zmax]');
% improvise outline of 3D box
plot3( [xmin,xmin], [ymin,ymin], [zmin,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymin,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmax], [ymin,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymin,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymin,ymin], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
xlabel('m_1');
ylabel('m_2');
zlabel('m_3');
% pretty crazy way to draw arrowheads!
% first arrowhead
rbar=m1;
plot3( [0, rbar(1)], [0, rbar(2)], [0, rbar(3)], 'k-', 'LineWidth', 3 );
tangent = rbar/sqrt(rbar'*rbar);
per1 = cross( tangent, [0, 0, 1]' );
per1 = per1/sqrt(per1'*per1);
per2 = cross( tangent, per1 );
per2 = per2/sqrt(per2'*per2);
L = 0.05;
v1 = rbar - L*tangent + 0.25*L*per1;
v2 = rbar - L*tangent - 0.25*L*per1;
v3 = rbar - L*tangent + 0.25*L*per2;
v4 = rbar - L*tangent - 0.25*L*per2;
plot3( [rbar(1), v1(1)], [rbar(2), v1(2)], [rbar(3), v1(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v2(1)], [rbar(2), v2(2)], [rbar(3), v2(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v3(1)], [rbar(2), v3(2)], [rbar(3), v3(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v4(1)], [rbar(2), v4(2)], [rbar(3), v4(3)], 'k-', 'LineWidth', 3 );
% second arrowhead
rbar=m2;
plot3( [0, rbar(1)], [0, rbar(2)], [0, rbar(3)], 'k-', 'LineWidth', 3 );
tangent = rbar/sqrt(rbar'*rbar);
per1 = cross( tangent, [0, 0, 1]' );
per1 = per1/sqrt(per1'*per1);
per2 = cross( tangent, per1 );
per2 = per2/sqrt(per2'*per2);
L = 0.05;
v1 = rbar - L*tangent + 0.25*L*per1;
v2 = rbar - L*tangent - 0.25*L*per1;
v3 = rbar - L*tangent + 0.25*L*per2;
v4 = rbar - L*tangent - 0.25*L*per2;
plot3( [rbar(1), v1(1)], [rbar(2), v1(2)], [rbar(3), v1(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v2(1)], [rbar(2), v2(2)], [rbar(3), v2(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v3(1)], [rbar(2), v3(2)], [rbar(3), v3(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v4(1)], [rbar(2), v4(2)], [rbar(3), v4(3)], 'k-', 'LineWidth', 3 );
% third arrowhead
rbar=m3;
plot3( [0, rbar(1)], [0, rbar(2)], [0, rbar(3)], 'k-', 'LineWidth', 3 );
tangent = rbar/sqrt(rbar'*rbar);
per1 = cross( tangent, [0, 0, 1]' );
per1 = per1/sqrt(per1'*per1);
per2 = cross( tangent, per1 );
per2 = per2/sqrt(per2'*per2);
L = 0.05;
v1 = rbar - L*tangent + 0.25*L*per1;
v2 = rbar - L*tangent - 0.25*L*per1;
v3 = rbar - L*tangent + 0.25*L*per2;
v4 = rbar - L*tangent - 0.25*L*per2;
plot3( [rbar(1), v1(1)], [rbar(2), v1(2)], [rbar(3), v1(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v2(1)], [rbar(2), v2(2)], [rbar(3), v2(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v3(1)], [rbar(2), v3(2)], [rbar(3), v3(3)], 'k-', 'LineWidth', 3 );
plot3( [rbar(1), v4(1)], [rbar(2), v4(2)], [rbar(3), v4(3)], 'k-', 'LineWidth', 3 );
% draw triangle connecting vector tips
plot3( [m1(1), m2(1)] , [m1(2), m2(2)], [m1(3), m2(3)], 'r:', 'LineWidth', 2 );
plot3( [m2(1), m3(1)] , [m2(2), m3(2)], [m2(3), m3(3)], 'r:', 'LineWidth', 2 );
plot3( [m3(1), m1(1)] , [m3(2), m1(2)], [m3(3), m1(3)], 'r:', 'LineWidth', 2 );
% set view
view(3);
fprintf('Caption: The three vectors (black) line in a plane (red), and so do not span the 3D space\n');
Caption: The three vectors (black) line in a plane (red), and so do not span the 3D space
% gdama09_03
% Householder rotation (over-determined problem)
% Note: No special effort is made to optimize the computations
disp('gdapy09_03');
gdapy09_03
% random matrix and data
N=10;
M=6;
G = random('uniform',-1.0,1.0,N,M);
mtrue = random('uniform',-1.0,1.0,M,1);
dobs = G*mtrue + random('Normal',0.0,0.1,N,1);
gda_draw(G,'caption G',dobs,' ', 'caption d');
% least squares olution for reference
GTG = G'*G;
GTd = G'*dobs;
mls = GTG\GTd;
els = dobs - G*mls;
Els = els'*els;
% setup for recursion
d = dobs;
T = zeros(N,N,M);
for i = (1:M)
a = sqrt(sum(G(i:N,i).^2));
a
if( G(i,i)<0.0 )
a=-a;
end
v = zeros(N,1);
v(i:N) = G(i:N,i);
v(i) = G(i,i)-a;
T(:,:,i) = eye(N)-(2.0*v*v')/(v'*v);
Ti=squeeze(T(:,:,i));
G = Ti*G;
d = Ti*d;
gda_draw(G,sprintf('caption G%d',i), ' ', d,sprintf('caption d%d',i) );
end
a = 1.6973
a = 1.9743
a = 1.9747
a = 1.7916
a = 0.6145
a = 1.3393
E = sum(d(M+1:N).^2);
fprintf("Total Error: least squares %f Householder %f\n", Els, E);
Total Error: least squares 0.014932 Householder 0.014932
fprintf(" \n");
% backsolve
m = zeros(M,1);
for i = (M:-1:1)
m(i) = d(i);
for j = (i+1:M)
m(i) = m(i) - G(i,j)*m(j);
end
m(i) = m(i)/G(i,i);
end
fprintf('m-ls m-House error\n');
m-ls m-House error
D = [mls,m,abs(mls-m)];
disp(D);
-0.3104 -0.3104 0.0000 -0.9644 -0.9644 0.0000 0.0776 0.0776 0.0000 -0.3331 -0.3331 0.0000 -0.5550 -0.5550 0.0000 -0.6146 -0.6146 0.0000
fprintf(" \n");
fprintf(" \n");
% gdama09_04
% transformation that simplifies L=m'*wm*m to L=mp'*mp
clear all;
fprintf('gdama09_04');
gdama09_04
% simple inverse problem
N=50;
M=N;
m = random('Normal',0,1,M,1); % a model chosen randomly
% prior information of flatness
epsi = 0.1;
D = epsi*toeplitz([1;-1;zeros(N-2,1)],[1,zeros(1,N-1)]);
Wm = D'*D;
% weighted length of m
L = m'*Wm*m;
% transformtion
Wmroot = sqrtm(Wm);
mp = Wmroot * m;
% unweigthed length of mp
Lp = mp'*mp;
fprintf('Weighted length of original %.4f\n',L);
Weighted length of original 0.9299
fprintf('Unweighted length of transformed %.4f\n',Lp);
Unweighted length of transformed 0.9299
gda_draw(D,'caption D',' ', Wmroot, 'caption Wm^1/2');
fprintf('Caption: The first difference operator D (left) is not\n');
Caption: The first difference operator D (left) is not
fprintf('equal to sqrt(Wm) (right), as the former is anti-symmetric\n');
equal to sqrt(Wm) (right), as the former is anti-symmetric
fprintf('and the latter symmetric\n');
and the latter symmetric
% gdama09_05
% solved Gm=d using Singular Value Decomposition
% for a G that averages adjacent model pararmaters
% supports Section 7.6 and Figure 7.4
clear all;
fprintf('gdama09_05\n');
gdama09_05
% model parameters m(z) where z is an auxially variable
M=20;
zmin = 0.0;
zmax = 5.0;
Lz = zmax-zmin;
Dz = (zmax-zmin)/(M-1);
z = zmin + Dz*[0:M-1]';
% tre model parameter is a sine wave
mtrue = sin( 2*pi*z/Lz );
% data kernel avarages L adjacent model parameters
% but wraps around from the end to the beginning of
% the model
g1 = [1, 1, 1.1, 1, 1]';
L=length(g1);
g=[ g1', zeros(1,M-L) ]';
N=M;
G=zeros(N,M);
for i=(1:M)
G(i,:) = circshift(g,i-1)';
end
G=fliplr(G);
% observed data is true data plus random noise
sigmad=0.01;
dobs = G*mtrue + random('Normal',0.0,sigmad,N,1);
% singular value decomposition of data kernel
[U, L, V] = svd(G,0);
lambda = diag(L);
% selection of non-zero singular valies
p = 16;
% throw out part of decomposition associated with zero singular values
lambdap = lambda(1:p);
Up=U(:,1:p);
Vp=V(:,1:p);
% natural solution
mest = Vp*((Up'*dobs)./lambdap);
% damped least squares solution, for comparison
e2=1e-2;
mestDML = G'*((G*G'+e2*eye(N,N))\dobs);
figure(2);
clf;
% plot singular values
set(gca,'LineWidth',2);
hold on;
axis( [1, M, 0 max(lambda) ] );
plot( [1:M]', lambda, 'k-', 'LineWidth', 2 );
hold on;
plot( [1:M]', lambda, 'ko', 'LineWidth', 2 );
xlabel('i');
ylabel('L_i');
fprintf('Caption: Singular values of the data kernel G\n');
Caption: Singular values of the data kernel G
% draw dobs = G m for various m's
gda_draw(G,' ',mtrue,'=',' ',dobs,' ',' ',' ',' ',mtrue,' ',mest,' ',mestDML);
fprintf('Caption: (Left) Graphical depiction of the equation Gm=d. (right)\n');
Caption: (Left) Graphical depiction of the equation Gm=d. (right)
fprintf('True solution, natural solution and damped minimum length solution.\n');
True solution, natural solution and damped minimum length solution.
% gdama09_06
% Singular Value Decomposition of two G's
% Both G's averages adjacent model pararmaters
% but the second is a smoother average
clear all;
fprintf('gdama09_06\n');
gdama09_06
% model parameters m(z) where z is an auxially variable
M=20;
mmin = 0.0;
mmax=5;
Dm = (mmax-mmin)/(M-1);
m = mmin + Dm*[0:M-1]';
% first G is a straight average, but wraps around the
% end of the (zmin,zmax) interval
g1 = [1, 1, 1, 1, 1]';
L=length(g1);
g=[ g1', zeros(1,M-L) ]';
N=M;
G1=zeros(N,M);
for i=(1:M)
G1(i,:) = circshift(g,i-1)';
end
G1=fliplr(G1);
% Singular Value Decomposition of G
[U, S1, V] = svd(G1,0);
s1 = zeros(M,1);
sp = diag(S1);
p = length(sp);
s1 = [ sp', zeros(1,M-p) ]';
% second G is a smoother average, and does not wrap
cmin=0;
cmax=10*(1/M);
c = cmin + (cmax-cmin)*[0:M-1]'/(M-1);
G2 = exp( -(c * [1:M]) ) - 0.9*exp( -(2*c * [1:M]).^2 );
% draw the two G's
gda_draw(G1,' ',G2);
fprintf('Caption: (Left) First data kernel. (Right) Second data kernel.\n');
Caption: (Left) First data kernel. (Right) Second data kernel.
figure();
clf;
% plot singular values
subplot(2,1,1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [1, M, 0 max(s1) ] );
plot( [1:M]', s1, 'k-', 'LineWidth', 2 );
hold on;
plot( [1:M]', s1, 'ko', 'LineWidth', 2 );
xlabel('i');
ylabel('S_i');
% singular value decomposition
[U, S1, V] = svd(G2,0);
s1 = zeros(M,1);
sp = diag(S1);
p = length(sp);
s2 = [ sp', zeros(1,M-p) ]';
% plot singular values
subplot(2,1,2);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [1, M, 0 max(s2) ] );
plot( [1:M]', s2, 'k-', 'LineWidth', 2 );
hold on;
plot( [1:M]', s2, 'ko', 'LineWidth', 2 );
xlabel('i');
ylabel('S_i');
fprintf('Caption: (Top) Singular-values of the first data kernel.\n');
Caption: (Top) Singular-values of the first data kernel.
fprintf('Caption: (Bottom) Singular-values of the second data kernel.\n');
Caption: (Bottom) Singular-values of the second data kernel.
% gdama09_07
% exercizing formula associated the the derivation of the SVD
clear all;
fprintf('gdama09_07\n');
gdama09_07
Etotal=0.0; % accumulated error
% exemplary sizes; requires N>M, M>p,
N=19;
M=7;
p=4;
% random matrix with p non-zero singular values
G = [ random('Normal',0,1,p,M); zeros(N-p,M)];
% S and its eigenvalues and eigenvectors. Flip the
% order to put the singular values in descending order
S=[zeros(N,N), G; G', zeros(M,M)];
[VS, LStrue] = eig(S);
lStrue = diag(LStrue);
lStrue = flipud(lStrue);
VS = fliplr(VS);
% eigenvalues and eigenvectors of GTG
% eigenvalues in ascending order
GTG = G'*G;
[V, L2V] = eig(GTG);
l2V = diag(L2V);
V = fliplr(V);
l2V = flipud(l2V);
lV = sqrt(abs(l2V));
Vp = V(:,1:p);
V0 = V(:,p+1:M);
Vz = zeros(M,N-M); % fills out so [Vp, Vp, Vz] is MxN
lVp = lV(1:p);
% eigenvalues and eigenvectors of GGT
% eigenvalues in ascending order
GGT = G*G';
[U, L2U] = eig(GGT);
l2U = diag(L2U);
U=fliplr(U);
l2U = flipud(l2U);
lU = sqrt(abs(l2U));
Up = U(:,1:p);
U0 = U(:,p+1:M);
Uz = U(:,M+1:N);
lUp = lU(1:p);
% check that GTG and GGT have the same non-zero eigenvalues
e = lVp-lUp;
El = max(abs(e(:)));
Etotal = Etotal + El;
fprintf('Maximum deviation of lVp from lUp: %.6e\n',El);
Maximum deviation of lVp from lUp: 1.998401e-15
% build eigenvectors of S from U and V
W = [Up, U0, sqrt(2)*Uz, -U0, -fliplr(Up); Vp, V0, Vz, V0, fliplr(Vp)]/sqrt(2);
% check that W is an orthogonal matrix saitifying WWT = WTW = 2I
% the reason for the sqrt(2) is that if the u's and v's are of
% unit lenth, then you need to rescale them when you build the w's.
WWT=W*W';
e = WWT-eye(M+N);
E1 = max(abs(e(:)));
WTW=W'*W;
e = WTW-eye(M+N);
E2 = max(abs(e(:)));
fprintf('Maximum deviation of WWT from 2I: %.6e\n',E1);
Maximum deviation of WWT from 2I: 4.440892e-16
fprintf('Maximum deviation of WTW from 2I: %.6e\n',E2);
Maximum deviation of WTW from 2I: 5.551115e-16
Etotal = Etotal + E1;
Etotal = Etotal + E2;
% WT*S*W should be the eigenvalues, so check that they are,
% but resort then to ensure consistent order. There is
% a subtle point here: there is no guarantee that a (u,v)
% pair generated from (GTG and GGT) corresponds to a positive
% singular value - it could be negative, in which case
% (-u,v) generates the positive singular value. So one needs
% to diddle with signs to get the singular values to get the
% list into decending order.
lSest = diag(W'*S*W);
lSpre = [lUp; zeros(N+M-2*p,1); -lUp];
for i = (1:p)
i2 = (N+M+1)-i;
if( lSest(i) < 0.0 )
lSest(i) = -lSest(i);
W(1:N,i) = -W(1:N,i);
lSest(i2) = -lSest(i2);
W(1:N,i2) = -W(1:N,i2);
end
end
% recompute eigenvalues and check in the right order
count = 0;
lSest = diag(W'*S*W);
tol=1e-6;
for i=(2:N+M)
e=lSest(i)-lSest(i-1);
if(e>tol)
count=count+1;
fprintf('%d %e\n',i,e );
end
end
fprintf('%d eigenvalues are out of order\n',count);
0 eigenvalues are out of order
Etotal = Etotal + count;
% recheck orthonormality
WWT=W*W';
e = WWT-eye(M+N);
E1 = max(abs(e(:)));
WTW=W'*W;
e = WTW-eye(M+N);
E2 = max(abs(e(:)));
fprintf('Maximum deviation of WWT from 2I: %.6e (2nd pass)\n',E1);
Maximum deviation of WWT from 2I: 4.440892e-16 (2nd pass)
fprintf('Maximum deviation of WTW from 2I: %.6e (2nd pass)\n',E2);
Maximum deviation of WTW from 2I: 5.551115e-16 (2nd pass)
Etotal = Etotal + E1;
Etotal = Etotal + E2;
e = lStrue - lSest;
ElS1 = max(abs(e));
fprintf('Maximum deviation of WT*S*W from eigenvalues: %.6e\n',ElS1);
Maximum deviation of WT*S*W from eigenvalues: 3.108624e-15
Etotal = Etotal + ElS1;
% [lUP, 0, -LUp] should be the eigenvalues, too,
% so check that they are, but resort then to ensure consistent order
lSpre = sort(lSpre,'descend');
e = lStrue - lSpre;
ElS2 = max(abs(e));
fprintf('Maximum deviation of (U,V) from eigenvalues: %.6e\n',ElS2);
Maximum deviation of (U,V) from eigenvalues: 1.776357e-15
Etotal = Etotal + ElS2;
% singular value decomposition
[U, L, V] = svd(G,0);
l = diag(L);
Vp = V(:,1:p);
Up = U(:,1:p);
lp = l(1:p);
% check singular vaues
e = lp - lSpre(1:p);
ElS3 = max(abs(e));
fprintf('Maximum deviation of SVD from eigenvalues: %.6e\n',ElS3);
Maximum deviation of SVD from eigenvalues: 1.776357e-15
Etotal = Etotal + ElS3;
% check reconstruction
Gest = U*L*V';
e = G(:)-Gest(:);
EG = max(abs(e))
EG = 8.8818e-16
fprintf('Maximum deviation of ULVT from G: %.6e\n',EG);
Maximum deviation of ULVT from G: 8.881784e-16
Etotal = Etotal + EG;
% check alternate reconstruction
Gest = Up*(Vp'.*lp);
e = G(:)-Gest(:);
EG = max(abs(e))
EG = 8.8818e-16
fprintf('Maximum deviation of ULVT from G: %.6e (method 2)\n',EG);
Maximum deviation of ULVT from G: 8.881784e-16 (method 2)
Etotal = Etotal + EG;
% warning: check will not work right if any non-zero singular values
% are degenerate
Eu = 0.0;
for i=(1:p)
uw = sqrt(2)*W(1:N,i);
usvd = U(:,i);
e = abs(1.0-abs(uw'*usvd));
Eu = Eu + e;
end
fprintf('Maximum deviation of U from Usvd: %.6e\n',Eu);
Maximum deviation of U from Usvd: 5.551115e-16
Etotal = Etotal + Eu;
% warning: check will not work right if any non-zero singular values
% are degenerate
Ev = 0.0;
for i=(1:p)
vw = sqrt(2)*W(N+1:N+M,i);
vsvd = V(:,i);
e = abs(1.0-abs(vw'*vsvd));
Ev = Ev + e;
end
fprintf('Maximum deviation of V from Vsvd: %.6e\n',Ev);
Maximum deviation of V from Vsvd: 1.221245e-15
Etotal = Etotal + Eu;
% if the total error is close to zero, then all tests
% succeeded and the derivation method really seems to
% reproduce the result of the SVD
fprintf('Overall Error %e\n',Etotal);
Overall Error 1.354472e-14
% gdama09_08
% graphical interpretation of Kuhn-Tucker theorem
clear all;
fprintf('gdama09_08\n');
gdama09_08
% (x,y) grid
N=51;
M=51;
Dx=0.01;
Dy=0.01;
x=Dx*[0:N-1]';
y=Dy*[0:M-1]';
% a long-wavelngth, 2D sinusoid
E=-sin(x)*sin(y)';
% minus gradient
dEdx = cos(x)*sin(y)';
dEdy = sin(x)*cos(y)';
% gradient list
gs = 0.1; % scale factor
gi = [40, 30, 40, 25, 40, 10, 10, 20, 10]';
gj = [40, 30, 25, 40, 10, 40, 10, 10, 20]';
Ng = length(gi);
% plot
figure(1);
clf;
set(gca,'LineWidth',3);
hold on;
set(gca,'FontSize',14);
colormap('jet');
axis ij;
axis( [x(1), x(N), y(1), y(M)] );
imagesc( [x(1), x(N)], [y(1), y(M)], E );
for k=(1:Ng)
i = gi(k);
j = gj(k);
plot( x(i), y(j), 'wo', 'LineWidth', 3 );
plot( [x(i),x(i)+gs*dEdx(i,j)]', [y(j),y(j)+gs*dEdy(i,j)]', 'w-', 'LineWidth', 3 );
end
% plot inequality constraint line
k1=7;
k2 = N-7+1;
plot( [x(k1), x(k2)]', [y(end), y(1)]', 'k-', 'LineWidth', 3);
% plot arrows on feasible side of equality
p1 = [x(k1), y(end)]';
p2 = [x(k2), y(1) ]';
t = p1-p2;
t = t/sqrt(t'*t);
n = [-t(2), t(1)]';
n = n/sqrt(n'*n);
sn = 0.03;
Nn = 10;
for k=(2:Nn-1)
a = (k-1)/(Nn-1);
p = a*p1 + (1-a)*p2;
plot( [p(1),p(1)+sn*n(1)]', [p(2),p(2)+sn*n(2)]', '-', 'LineWidth',2, 'Color', [0.5,0.5,0.5] );
end
% tabulate the deviation between n and -gradE/|gradE| along
% the constraint line
Nn = 101;
e = zeros(Nn,1);
ie = zeros(Nn,1);
je = zeros(Nn,1);
for k=(1:Nn)
a = (k-1)/(Nn-1);
p = a*p1 + (1-a)*p2;
i = floor(p(1)/Dx)+1;
j = floor(p(2)/Dy)+1;
g = -[dEdx(i,j), dEdy(i,j)]';
g = g/sqrt(g'*g);
e(k) = (n(1)-g(1))^2 + (n(2)-g(2))^2;
ie(k) = i;
je(k) = j;
end
% find and plot point where gradE and n are anti-parallel
[emin, k] = min(e);
i = ie(k);
j = je(k);
plot( x(i), y(j), 'ko', 'LineWidth', 4 );
xlabel('m1');
ylabel('m2');
colorbar();
fprintf('Caption: Graphical interpretation of the Kuhn-Tucker theorem, showing\n');
Caption: Graphical interpretation of the Kuhn-Tucker theorem, showing
fprintf('error surface (colors) and its gradient (while lines with dot on uphill end),\n');
error surface (colors) and its gradient (while lines with dot on uphill end),
fprintf('constraint surface (black line) and its feasible-pointing normals (grey bars), and fessible point\n');
constraint surface (black line) and its feasible-pointing normals (grey bars), and fessible point
fprintf('with minimum error.\n');
with minimum error.
% gdama09_09
% non-negative least squares applied to a simple straight line problem
clear all;
fprintf('gdama09_09\n');
gdama09_09
% z-axis
N=11;
zmin = 0.0;
zmax=10.0;
Dz = (zmax-zmin)/(N-1);
z = zmin + Dz*[0:N-1]';
% data kernel
M = 2;
G = [ones(N,1), z];
% synthetic data
mtrue = [-2, 1]';
dtrue = G*mtrue;
sd=2.0;
dobs = dtrue + random('Normal',0,sd,N,1);
% ordinary least square solution for comparison
mestls = (G'*G)\(G'*dobs);
dprels = G*mestls;
e = (dobs-dprels);
Els=e'*e;
% non-negative least squares
mest = lsqnonneg(G,dobs);
dpre = G*mest;
E = (dobs-dpre)'*(dobs-dpre);
figure();
clf;
hold on;
set(gca,'LineWidth',2);
set(gca,'FontSize',14);
axis( [zmin, zmax, -10, 20]);
plot(z,dtrue,'k-','LineWidth',4);
plot(z,dobs,'b*','LineWidth',2);
plot(z,dprels,'r-','LineWidth',2);
plot(z,dpre,'g-','LineWidth',2);
xlabel('z');
ylabel('d');
fprintf('Caption: Data (circles) fit by strsight lines, true (black), least squares\n');
Caption: Data (circles) fit by strsight lines, true (black), least squares
fprintf('(red), non-negative least squares (green)\n');
(red), non-negative least squares (green)
fprintf('\n');
fprintf('mtrue %.2f %.2f Error %.2f\n', mtrue(1), mtrue(2), 0.0);
mtrue -2.00 1.00 Error 0.00
fprintf('mls %.2f %.2f Error %.2f\n', mestls(1), mestls(2), Els);
mls -4.09 1.36 Error 14.35
fprintf('mnnls %.2f %.2f Error %.2f\n', mest(1), mest(2), E);
mnnls 0.00 0.78 Error 66.87
% gdama09_10
% illustration of natural and non-negative least squares
% for gravity problem
clear all;
fprintf('gdama09_10\n');
gdama09_10
% model is two-dimensional grid of pixels, each of which represents mass
% coorinate system, x is down, y is right, origin in upper-left
Nx = 20;
Ny = 20;
M=Nx*Ny;
% model x coordinate
xmin = 0.0;
xmax = 1.0;
Lx = xmax-xmin;
Dx = (xmax-xmin)/(Nx-1);
x = xmin + Dx*[0:Nx-1]';
% model y coordinate
ymin = 0.0;
ymax = 1.0;
Ly = ymax-ymin;
Dy = (ymax-ymin)/(Ny-1);
y = ymin + Dy*[0:Ny-1]';
% gravity data are on top surface, extending past edges of model by several
% model widths
N=30*Ny;
xo = -0.5; % observations well above surface of model
yomin=ymin-2*Ly;
yomax=ymax+2*Ly;
Lyo = (yomax-yomin);
Dyo = (yomax-yomin)/(N-1);
yo = yomin + Dyo*[0:N-1]';
% data kernel for vertical component of gravity
g = 1; % gravitational constant (not a relistic value)
G=zeros(N,M);
jofq=zeros(Nx,1); % backpointer to x(j) value of model parameters
kofq=zeros(Ny,1); % backpointer to y(k) value of model parameters
for i = (1:N) % loop over observation points
% loop over masses
q = 0; % counts model parameters
for j = (1:Nx) % x
for k = (1:Ny) % y
q = q+1;
jofq(q)=j; kofq(q)=k;
D2 = (x(j)-xo)^2 + (y(k)-yo(i))^2; % distance squared
ct = (x(j)-xo) / sqrt(D2); % cosine of vertical angle
G(i,q) = ct*g/D2;
end
end
end
% true model parameters; just something cooked up to look interesting
mxytrue = ([1:Nx]-Nx/2)'*ones(1,Ny)/Nx * ones(Nx,1)*([1:Ny]-Ny/2)/Ny;
mxytrue = mxytrue - min(min(mxytrue));
mxytrue = (mxytrue / max(max(mxytrue))) + random('uniform',0,0.1,Nx,Ny);
% convert to vector
mtrue=zeros(M,1);
for q = (1:M)
mtrue(q) = mxytrue( jofq(q), kofq(q) );
end
% convert back to image as a check
% mxytrue2=zeros(Nx,Ny);
% for q = [1:M]
% mxytrue2( jofq(q), kofq(q) ) = mtrue(q);
% end
% gda_draw(' ',mxytrue,' ',mxytrue2);
% create and plot synthetic data
sigmad=5;
dobs = G*mtrue + random('Normal',0.0,sigmad,N,1);
% calculate and plot singular values
[U, L, V] = svd(G,0);
lambda = diag(L);
p=15;
lambdap=lambda(1:p);
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [1, length(lambda), 0 max(lambda) ] );
plot( [1:length(lambda)]', lambda, 'k-', 'LineWidth', 2 );
plot( [1:length(lambda)]', lambda, 'ko', 'LineWidth', 2 );
plot( [1:p]', lambdap, 'ro', 'LineWidth', 2 );
xlabel('i');
ylabel('L_i');
fprintf('Caption: Singular values\n');
Caption: Singular values
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [yomin, yomax, 0, max(dobs) ] );
plot( yo, dobs, 'k-', 'LineWidth', 3 );
xlabel('y');
ylabel('data d(y)');
% natural solution
Up=U(:,1:p);
Vp=V(:,1:p);
lambdap=lambda(1:p);
mestN = Vp*((Up'*dobs)./lambdap);
% convert to image
mxyestN=zeros(Nx,Ny);
for q = (1:M)
mxyestN( jofq(q), kofq(q) ) = mestN(q);
end
% natural solution, smaller p
p=4;
Up=U(:,1:p);
Vp=V(:,1:p);
lambdap=lambda(1:p);
mestN2 = Vp*((Up'*dobs)./lambdap);
% convert to image
mxyestN2=zeros(Nx,Ny);
for q = (1:M)
mxyestN2( jofq(q), kofq(q) ) = mestN2(q);
end
% predicted data
dpreN = G*mestN;
plot( yo, dpreN, 'ro', 'LineWidth', 2 );
eN = dobs-dpreN;
EN=eN'*eN;
% non-negative splution
mestNN = lsqnonneg(G,dobs);
% convert to image
mxyestNN=zeros(Nx,Ny);
for q = (1:M)
mxyestNN( jofq(q), kofq(q) ) = mestNN(q);
end
% predicted data
dpreNN = G*mestNN;
plot( yo, dpreNN, 'g.', 'LineWidth', 2 );
eNN = dobs-dpreNN;
ENN=eNN'*eNN;
fprintf('Caption: Observed data (black), data predicted from natural solution (red),\n');
Caption: Observed data (black), data predicted from natural solution (red),
fprintf('data predicted from non-negative solution (green)\n');
data predicted from non-negative solution (green)
gda_draw(' ',mxytrue,' ',mxyestN2,' ',mxyestN,' ',mxyestNN);
fprintf('Caption: (Left) True solution, (2nd) Natural solution with small p,\n');
Caption: (Left) True solution, (2nd) Natural solution with small p,
fprintf('(3rd) Natural solution with large p, (right) Non-negative solution\n');
(3rd) Natural solution with large p, (right) Non-negative solution
fprintf('Total Error: Natural %f Non-Negative %f\n',EN,ENN);
Total Error: Natural 13810.460981 Non-Negative 13859.392103
% gdama09_11
% depiction of feasible areas, for 4 cases
clear all;
fprintf('gdama09_11\n');
gdama09_11
m1min=0;
m1max=2;
Nm1 = 41;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = m1min + Dm1*[0:Nm1-1]';
m2min=0;
m2max=2;
Nm2 = 41;
Dm2 =(m2max-m2min)/(Nm2-1);
m2 = m2min + Dm2*[0:Nm2-1]';
X = m1*ones(Nm2,1)';
Y = ones(Nm1,1)*m2';
% three inequality constraints
% m1-m2>=0 [1 -1] [0]
% 0.5*m1+m2>=1 [0.5 1] [1]
% m1>=0.2 [1 0] [0.2]
% in the last constraint, the 0.2 can be diddled
% right now this constraint does not affect the
% feasible region, but changing 0.2 to 1.2 changes
% the region considerably
H = [ [1, 0.5, 1]', [-1, 1, 0]' ];
h = [0, 1, 0.2]';
p=1;
C1 = zeros(Nm1,Nm2);
for i=(1:Nm1)
for j=(1:Nm2)
m = [ m1(i), m2(j) ]';
C1(i,j)=H(p,:)*m-h(p);
end
end
C1a = (C1>=0)+0.001;
p=2;
C2 = zeros(Nm1,Nm2);
for i=(1:Nm1)
for j=(1:Nm2)
m = [ m1(i), m2(j) ]';
C2(i,j)=H(p,:)*m-h(p);
end
end
C2a = (C2>=0)+0.001;
p=3;
C3 = zeros(Nm1,Nm2);
for i=(1:Nm1)
for j=(1:Nm2)
m = [ m1(i), m2(j) ]';
C3(i,j)=H(p,:)*m-h(p);
end
end
C3a = (C3>=0)+0.001;
CT = ((C1>=0)&(C2>=0)&(C3>=0))+0.001;
Gp = [H, h]';
dp = [zeros(1,length(H(1,:))), 1]';
mp = lsqnonneg(Gp,dp);
ep = dp - Gp*mp;
m = -ep(1:end-1)/ep(end);
im1est = floor((m(1)-m1min)/Dm1)+1;
jm2est = floor((m(2)-m2min)/Dm2)+1;
CT(im1est,jm2est)=0.5;
gda_draw(C1a,'caption (A)',' ',C2a,'caption (B)',' ',C3a,'caption (C)',' ',CT,'caption (D)');
fprintf('Caption: (Left three frames) Three inequality constraints, esch dividing\n');
Caption: (Left three frames) Three inequality constraints, esch dividing
fprintf('the (m2,m2) plane into feasible (brown) and infeasible (blue( region. (Right)\n');
the (m2,m2) plane into feasible (brown) and infeasible (blue( region. (Right)
fprintf('Union of the three inequality constraints, with green dot the feasible point\n');
Union of the three inequality constraints, with green dot the feasible point
fprintf('that minimizes |m|^2\n');
that minimizes |m|^2
fprintf('\n');
e=H*m-h;
for i=(1:3)
fprintf('Constraint %d: Hm-h = %f\n', i, e(i) );
end
Constraint 1: Hm-h = 0.000000 Constraint 2: Hm-h = -0.000000 Constraint 3: Hm-h = 0.466667
% gdama09_12
% least squares with inequality constraints
% exemplary stright line problem
clear all;
fprintf('gdama09_12\n');
gdama09_12
% three inequality constraints
% m1-m2>=0 [1 -1] [0]
% 0.5*m1+m2>=1 [0.5 1] [1]
% m1>=0.2 [1 0] [0.2]
% in the last constraint, the 0.2 can be diddled
% right now this constraint does not affect the
% feasible region, but changing 0.2 to 1.2 cnages
% the region considerably
H = [ [1, 0.5, 1]', [-1, 1, 0]' ];
h = [0, 1, 0.2]';
% straight line problem
% in feasible region: mtrue = [1, 0.8]';
% out of feasbile region mtrue = [1, 1.3]';
mtrue = [1.0, 0.2]';
N=11;
z = [0:N-1]'/(N-1);
G = [ones(N,1), z];
dtrue = G*mtrue;
sd=0.05;
dobs = dtrue + random('Normal',0,sd,N,1);
% solve by ordinary least squares, just for check
mestls = (G'*G)\(G'*dobs);
dprels = G*mestls;
fprintf('ordinary least squares:\n');
ordinary least squares:
fprintf(' mtrue %f %f\n', mtrue(1), mtrue(2) );
mtrue 1.000000 0.200000
fprintf(' mest %f %f\n', mestls(1), mestls(2) );
mest 0.989730 0.183593
fprintf('\n');
% (m1, m2) grid for plotting purposes
m1min=0;
m1max=2;
Nm1 = 41;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = m1min + Dm1*[0:Nm1-1]';
m2min=0;
m2max=2;
Nm2 = 41;
Dm2 =(m2max-m2min)/(Nm2-1);
m2 = m2min + Dm2*[0:Nm2-1]';
X = m1*ones(Nm2,1)';
Y = ones(Nm1,1)*m2';
% Error surface, for plotting purposes
E = zeros(Nm1,Nm2);
for i=(1:Nm1)
for j=(1:Nm2)
m = [ m1(i), m2(j) ]';
e = dobs - G*m;
E(i,j)=e'*e;
end
end
Emin = min(min(E));
Emax = max(max(E));
% evaluate the constraints, for plotting purposes
% constraint 1
p=1;
C1 = zeros(Nm1,Nm2);
for i=(1:Nm1)
for j=(1:Nm2)
m = [ m1(i), m2(j) ]';
C1(i,j)=H(p,:)*m-h(p);
end
end
C1a = (C1>=0)+0.001;
% constraint 2
p=2;
C2 = zeros(Nm1,Nm2);
for i=(1:Nm1)
for j=(1:Nm2)
m = [ m1(i), m2(j) ]';
C2(i,j)=H(p,:)*m-h(p);
end
end
C2a = (C2>=0)+0.001;
% constraint 3
p=3;
C3 = zeros(Nm1,Nm2);
for i=(1:Nm1)
for j=(1:Nm2)
m = [ m1(i), m2(j) ]';
C3(i,j)=H(p,:)*m-h(p);
end
end
C3a = (C3>=0)+0.001;
% union of constraints
CT = ((C1>=0)&(C2>=0)&(C3>=0))+0.001;
% plot true solution by making dot in CT
im1est = floor((mtrue(1)-m1min)/Dm1)+1;
jm2est = floor((mtrue(2)-m2min)/Dm2)+1;
CT(im1est,jm2est)=0.7;
E(im1est,jm2est)=0.7*(Emax-Emin);
% SVD
% note that G must be overdetermined so V=Vp
% economy size (flag 0) SVD so U is Up
[Up, Lp, Vp] = svd(G,0);
lambda = diag(Lp);
rlambda = 1./lambda;
Lpi = diag(rlambda);
% transformation 1
Hp = -H*Vp*Lpi;
hp = h + Hp*Up'*dobs;
% transformation 2
Gp = [Hp, hp]';
dp = [zeros(1,length(Hp(1,:))), 1]';
mpp = lsqnonneg(Gp,dp);
ep = dp - Gp*mpp;
mp = -ep(1:end-1)/ep(end);
% take mp back to m
mest = Vp*Lpi*(Up'*dobs-mp);
dpre = G*mest;
% plot solution by making dot in CT
im1est = floor((mest(1)-m1min)/Dm1)+1;
jm2est = floor((mest(2)-m2min)/Dm2)+1;
CT(im1est,jm2est)=0.5;
E(im1est,jm2est)=0.5*(Emax-Emin);
gda_draw(' ',CT,'caption (A)',' ',E,'caption (B)');
fprintf('(left) Feasible region (brown) with error surface (colors), unconstrained\n');
(left) Feasible region (brown) with error surface (colors), unconstrained
fprintf('solution (orange dot) and constrained solution (green dot). (Right) Same as\n');
solution (orange dot) and constrained solution (green dot). (Right) Same as
fprintf('left, but with feasible region not shown\n');
left, but with feasible region not shown
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, 1, 0, 2 ] );
plot( z, dobs, 'ko', 'LineWidth', 3 );
plot( z, dtrue, 'k-', 'LineWidth', 3);
plot( z, dprels, 'r-', 'LineWidth', 3);
plot( z, dpre, 'g-', 'LineWidth', 3);
fprintf('Caption: Straight line fit, showing true data (black line), observed data\n');
Caption: Straight line fit, showing true data (black line), observed data
fprintf('(circles), unconstrained least squares fit (red), and constraind fit (green).\n');
(circles), unconstrained least squares fit (red), and constraind fit (green).
fprintf('least squares with inequality constraint:\n');
least squares with inequality constraint:
fprintf(' mtrue %f %f\n', mtrue(1), mtrue(2) );
mtrue 1.000000 0.200000
fprintf(' mest %f %f\n', m(1), mestls(2) );
mest 2.000000 0.183593
fprintf('\n');