Matlab code for this type of factorization.A has a SR decomposition A = SR , where S ∈ R^ 2n ×2n is a symplectic matrix, i.e. S ^TJS = J

Question

Farooq Aamir 2018-2-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/381885-matlab-code-for-this-type-of-factorization-a-has-a-sr-decomposition-a-sr-where-s-r-2n-x2n-i

回答： Farooq Aamir 2018-3-20

J=[0 -I;I 0] where I∈R^ n ×n means identity matrix. R=[R11 R12;R21 R22]∈ R^ 2n ×2n , is constituted by upper triangular matrices R11 , R12 , R22 and strictly upper triangular matrix R21. diag (R11 ) =|diag (R22 )| and diag (R12 ) = 0 , then the SR decomposition is unique.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Christine Tobler 2018-2-12

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/381885-matlab-code-for-this-type-of-factorization-a-has-a-sr-decomposition-a-sr-where-s-r-2n-x2n-i#answer_304780

I'm not very acquainted with the SR decomposition (all I know about it I found just now by googling). I would suggest to take a look at what seems to be the original paper, "Matrix factorizations for symplectic QR-like methods" by Bunse-Gerstner. This seems to suggest an algorithm using other factorizations (e.g. QR) as building blocks.

I'm also tagging this Control, because it seems that this decomposition has applications in total control.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Christine Tobler 2018-2-12

Also, take a look at the acknowledgements here: "On the sensitivity of the SR decomposition", Xiao-Wen Chang. They mention some MATLAB code for computing the SR decomposition being shared between researchers. Contacting one of them directly might be your best bet (although the paper is from 1998, so the code may not be easy to find).

请先登录，再进行评论。

Answer 2

Farooq Aamir 2018-3-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/381885-matlab-code-for-this-type-of-factorization-a-has-a-sr-decomposition-a-sr-where-s-r-2n-x2n-i#answer_310813

在 MATLAB Online 中打开

if true
  % code
end
function [c, v] = optsymhouse1(a)
twon = length(a); n = twon/2;
J = [zeros(n), eye(n); -eye(n), zeros(n)];
p = sign(a(1))*norm(a); aux = a(1)-p ;
if aux == 0
c = 0; v = zeros(twon, 1); %T = eye(twon);
elseif a(n + 1) == 0
display('division by zero');
return
else
v =a/aux;
c =aux^2/(p * a(n + 1));
v(1) = 1;
%T = (eye(twon) + c * v * v' * J);
end
function [c, v] =optsymhouse2(u)
twon = length(u); n = twon/2;
J = [zeros(n), eye(n); -eye(n), zeros(n)];
if n == 1
v = zeros(twon, 1); c = 0; %T = eye(twon);
else
I = [2 : n, n + 2 : twon]; e = norm(u(I));
if e == 0
v = zeros(twon, 1); c = 0; %T = eye(twon);
else
v1 = u(n + 1);
if v1 == 0
display('division by zero')
return
else
v = -u/e; v(1) = 1; v(n + 1) = 0; c = e/v1;
T = (eye(twon) + c * v * v' * J);
end
end
end
function[S,R] = SROSH(A)
[den, dep] = size(A);
n = den/2; p = dep/2;
v = zeros(den, 1);
JJ = [zeros(n) eye(n); -eye(n) zeros(n)];
S = eye(den);
for k = 1 : p
J= [zeros(n - k + 1), eye(n - k + 1); -eye(n - k + 1), zeros(n - k + 1)];
% Computing T 2k?1:
a = A([k : n, n + k : den],[k]);
[ c3,v3] = optsymhouse1(a);
T = eye(den-2*k+2) + c3*(v3*v3')*J ; % T not formed explicitly.
% Updating A:
A([k : n, n + k : den], [k : p, p + k : dep]) = T * A([k : n, n + k : den], [k : p, p + k :dep]);
% Computing T 2k?1 :
v(k : n) = v3( 1 : n - k + 1);
v(k + n : den) = v3(n -k + 2 : 2 * (n - k + 1));
TJ = eye(den) - c3*v*v'*JJ; % T J not formed explicitly;
% Updating S:
S = S * TJ;
% Computing T 2k:
u = A([k : n, n + k : den], p + k);
[ c3,v3] = optsymhouse2(u);
T = eye(den -2*k + 2) + c3 * (v3*v3')*J ;  % T not formed explicitly.
% Updating A[:
A([k : n, n + k : den], [k : p, p + k : dep]) = T *A([k : n, n + k : den], [k : p, p + k :dep]);
% Computing T 2k:
v(k  : n) = v3(1 : n - k + 1);
v(k  + n : den) = v3(n  - k + 2 : 2 * (n - k + 1));
TJ1 = (eye(den) - c3*(v*v')*JJ); % T J not formed explicitly;
% Updating S:
S = S * TJ1;
end
R = A;
%JJ'*S'*JJ*S = eye(2n);