gusucode.com > qit_matlab_0.10.0工具箱源码程序 > qit/utils/superop_fp.m
function [A, spectrum] = superop_fp(L, rtol) % SUPEROP_FP Fixed point states of a Liouvillian superoperator. % [A, spectrum] = superop_fp(L, tol) % % Finds the intersection of the kernel of the Liouvillian L % with the set of valid state operators, thus giving the set of % fixed point states for the quantum channel represented by the % master equation % % \dot{\rho} = inv_vec(L * vec(\rho)). % % Let size(L) == [D, D] (and d = sqrt(D) be the dimension of the Hilbert space). % % Returns the D * n matrix A, which contains as its columns a set % of n vectorized orthogonal Hermitian matrices (with respect to % the Hilbert-Schmidt inner product) which "span" the set of FP % states in the following sense: % % vec(\rho) = A * c, where c \in \R^n and c_1 = 1. % % A(:,1) is the shortest vector in the Hermitian kernel of L that % has trace 1, the other columns of A are traceless and normalized. % Hence, A defines an (n-1)-dimensional hyperplane. % % A valid state operator also has to fulfill \rho \ge 0. These % operators form a convex set in the Hermitian trace-1 hyperplane % defined by A. Currently this function does nothing to enforce % positivity, it is up to the user to choose the coefficients a_k % such that this condition is satisfied. % % Singular values of L less than or equal to the tolerance tol are % treated as zero. % If L has Lindblad form, if L(rho) = \lambda * rho, % we also have L(rho') = conj(\lambda) * rho' % Hence the non-real eigenvalues of L come in conjugate pairs. % % Especially if rho \in Ker L, also rho' \in Ker L. % Ville Bergholm 2011-2014 if nargin < 2 rtol = eps(class(L)); end % Hilbert space dimension d = sqrt(size(L, 2)); % Get the kernel of L, restricted to vec-mapped Hermitian matrices. % columns of A: complex orthogonal vectors A_i spanning the kernel (with real coefficients) [A, spectrum] = nullspace_hermitian(L, rtol); % Extract the trace-1 core of A and orthonormalize the rest. % We want A_i to be orthonormal wrt. the H-S inner product % <X, Y> := trace(X' * Y) = vec(X)' * vec(Y) % With this inner product, trace(A) = <eye(d), A>. temp = vec(speye(d)); a = (temp' * A)'; % Construct the shortest linear combination of A_i which has tr = 1. % This is simply (1/d) * eye(d) _iff_ it belongs to span(A_i). core = A * (a / norm(a)^2); % trace-1 % Remove the component parallel to core from all A_i. temp = core / norm(core); % normalized A = A -temp * (temp' * A); % Re-orthonormalize the vectors, add the core A = [core, orthonormalize(A, 1e-6)]; % FIXME tolerance N = size(A, 2); for k = 1:N temp = inv_vec(A(:, k)); A(:, k) = vec(0.5 * (temp + temp')); % re-Hermitize to fix numerical errors end % TODO intersect with positive ops %probe(A); %[u,t] = schur(L); %unnormality = norm(t, 'fro')^2 -norm(diag(t))^2 %[us, ts] = ordschur(u, t, 'udi') %E = ordeig(t); % TODO eigendecomposition, find orthogonal complement to span(v) = ker(v'). % these are the states which do not belong to eigenspaces and % should show transient polynomial behavior end function probe(A) % Probe the positivity of the FP set N = size(A, 2); % upper limit for amplitudes of orthogonal non-core state components (from purity): % necessary but not sufficient C = A(:, 1); a_max = sqrt(1 - C'*C); x = 1.1 * a_max; a = linspace(-x, x, 57); b = linspace(-x, x, 52); neg = []; for j=1:length(b) for k=1:length(a) r = C +a(k)*A(:, 2) +b(j)*A(:, 3); neg(j, k) = min(eig(inv_vec(r))); end end figure(); mesh(a, b, neg); xlabel('A_1') ylabel('A_2') zlabel('min eigenvalue') end function plot_eigenvalues(s) % Plots the eigenvalues of L on the complex plane. % propagate the eigenvalues a bit to see their behavior t = linspace(0, 1, 10); %for k = 1:length(t) % temp = expm(t(k) * L); % s_expL(:, k) = sort(eig(temp)); %end s_expL = exp(s * t); % FIXME assuming the normality/diagonalizability of L... if isreal(s_expL) s_expL = complex(s_expL); end figure(); % plot the unit circle phi = linspace(0, 2*pi, 51); plot(exp(1i * phi), '--'); %axis([-2 2 -2 2]); hold on plot(s, 'rx') plot(s_expL.', 'o-') legend('Unit circle', '\sigma(L)', '\sigma(exp(Lt))') title('Eigenvalue plot') axis equal grid on end