Skip to content
Snippets Groups Projects
Commit 43d56ae0 authored by Dan Waggoner's avatar Dan Waggoner
Browse files

modified and added some msv matlab files

parent b64f9fe0
No related branches found
No related tags found
No related merge requests found
...@@ -17,8 +17,10 @@ function [F1, F2, G1, G2, V, err] = fwz_msv_msre(P, A, B, Psi, s, x, max_count, ...@@ -17,8 +17,10 @@ function [F1, F2, G1, G2, V, err] = fwz_msv_msre(P, A, B, Psi, s, x, max_count,
% %
% eta(t) = F2{s(t)}*x(t-1) + G2{s(t)}*epsilont(t) % eta(t) = F2{s(t)}*x(t-1) + G2{s(t)}*epsilont(t)
% %
% err is set to 0 upon sucess. A positive value is the number of % A positive value of err is the number of iterations needed to obtain
% iterations before the method was terminated. % convergence and indicates success. A negitive value of err is the number of
% iterations before the method terminated without convergence and indicates
% failure.
% %
h=size(P,1); h=size(P,1);
...@@ -38,20 +40,9 @@ if nargin <= 5 ...@@ -38,20 +40,9 @@ if nargin <= 5
end end
f=zeros(h*s*(n-s),1); f=zeros(h*s*(n-s),1);
Is=eye(s); Is=eye(s);
I=eye(n-s); I=eye(n-s);
% The following code would be used if Pi were passed instead of the
% assumption that Pi' = [zeros(s,n-s) eye(s)].
% for i=1:h
% [Q,R] = qr(Pi{i});
% W=[zeros(n-s,s) I; inv(R(1:s,1:s)); zeros(s,n-s)]*Q';
% A{i}=W*A{i};
% B{i}=W*B{i};
% Psi{i}=W*Psi{i};
% end
U=cell(h,1); U=cell(h,1);
for i=1:h for i=1:h
U{i}=inv(A{i}); U{i}=inv(A{i});
...@@ -109,9 +100,9 @@ while cont ...@@ -109,9 +100,9 @@ while cont
end end
if (norm(f) < tol) if (norm(f) < tol)
err=0;
else
err=count; err=count;
else
err=-count;
end end
F1=cell(h,1); F1=cell(h,1);
......
function [A_new, B_new, Psi_new] = fwz_pi2eye(A, B, Psi, Pi)
% [A, B, Psi] = fwz_pi2eye(A, B, Psi, Pi)
% Converts the representation
%
% A(s(t) x(t) = B(s(t) x(t-1) + Psi(s(t)) epsilon(t) + Pi{s(t)} eta(t)
%
% to a form compatible with fwz_msv_msre(). This is the form
%
% A(s(t) x(t) = B(s(t) x(t-1) + Psi(s(t)) epsilon(t) + Pi eta(t)
%
% where Pi' = [zeros(s,n-s) eye(s)]
%
for i=1:h
[Q,R] = qr(Pi{i});
W=[zeros(n-s,s) I; inv(R(1:s,1:s)); zeros(s,n-s)]*Q';
A_new{i}=W*A{i};
B_new{i}=W*B{i};
Psi_new{i}=W*Psi{i};
end
\ No newline at end of file
function diff = verify_solution(P,A,B,Psi,s,F1,F2,G1,G2,V)
%
% Verifies that
%
% x(t) = V(s(t))(F1(s(t))x(t-1) + G1(s(t))epsilon(t)
% eta(t) = F2(s(t))x(t-1) + G2(s(t))epsilon(t)
%
% is a solution of
%
% A(s(t)) x(t) = B(s(t)) x(t-1) + Psi(s(t) epsilon(t) + Pi eta(t)
%
% where Pi' = [zeros(s,n-s) eye(s)] by verifying that
%
% [A(j)*V(j) Pi] [ F1(j)
% F2(j) ] = B(j)
%
% [A(j)*V(j) Pi] [ G1(j)
% G2(j) ] = Psi(j)
%
% (P(i,1)*F2(1) + ... + P(i,h)*F2(h))*V(i) = 0
%
h=size(P,1);
n=size(A{1},1);
r=size(Psi{1},2);
diff=0;
Pi=[zeros(n-s,s); eye(s)];
for j=1:h
X=inv([A{j}*V{j} Pi]);
tmp=norm(X*B{j} - [F1{j}; F2{j}]);
if tmp > diff
diff=tmp;
end
tmp=norm(X*Psi{j} - [G1{j}; G2{j}]);
if (tmp > diff)
diff=tmp;
end
end
for i=1:h
X=zeros(s,n);
for j=1:h
X=X+P(i,j)*F2{j};
end
tmp=norm(X*V{i});
if tmp > diff
diff=tmp;
end;
end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment