-
Notifications
You must be signed in to change notification settings - Fork 0
/
prepare_onesite_truncate.m
164 lines (154 loc) · 6.45 KB
/
prepare_onesite_truncate.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
function [B, U, para, results, sv, vNE] = prepare_onesite_truncate(A,para,sitej,results)
% Truncate Bond dimension via SV of A
%
% Changed:
% FS 20/10/2014: - removed explicit direction, using para.sweepto now!
% FS 28/10/2014: - implemented sweep to right for TDVP
% - Added keeping 1 more dim if smallest SV too large,
% limited by para.tdvp.maxBondDim in 'r' sweep
% - major restructuring using truncation function
% FS 09/04/2015: - in 'r' sweep: apply to site L to normalize MPS. No
% truncation or para.D saving!
[D1, D2, d] = size(A);
D_old = para.D;
err = 0;
switch para.sweepto
case 'r'
if para.parity=='n'
%% normal Parity
A = permute(A, [3, 1, 2]);
A = reshape(A, [d * D1, D2]); % reshapes to (a1 d),a2
[B, S, U] = svd2(A); % Anew with orthonormal columns
if norm(diag(S)) ~= 1
S = S./norm(diag(S)); % normalisation necessary for thermal steps
end
sv = diag(S);
if para.tdvp.truncateExpandBonds && sitej ~= para.L % only used in TDVP. For VMPS add: && sitej <= para.trustsite(end)
if sv(end) > para.svmintol
%% Expand A dims
if sv(end) > para.svmaxtol
dimincratio = log10(sv(end)/para.svmaxtol)/2; % aggressive expansion if min(SV) > SVmaxtol.
else
dimincratio = log10(sv(end)/para.svmintol)/(abs(log10(sv(end)))-1); % new scheme, also weighted by current lowest exponent, less aggressive than above
end
adddim = ceil(D2*dimincratio);
adddim = min(adddim, para.tdvp.maxBondDim(end) - para.D(sitej));
A = cat(2,A,zeros(d * D1,adddim));
[B, S, U] = svd2(A); % better since B contains no zeros but orthonormal vectors!
if norm(diag(S)) ~= 1
S = S./norm(diag(S)); % normalisation necessary for thermal steps
end
U = U(:,1:D2); % form old Bond Dim in columns == U=cat(2,U,addmat);
else
%% Truncate A dims
[B, S, U, err] = truncate(B,S,U);
end
end
DB = size(S, 1); % new a2 dimension of B
if sitej ~= para.L
para.D(sitej)=DB;
end
B = reshape(B, [d, D1, DB]);
B = permute(B, [2, 3, 1]);
U = S * U;
else
error('Have not implemented different parity yet!');
end
resultsSite = sitej;
case 'l'
if para.parity=='n'
A = permute(A, [1, 3, 2]);
A = reshape(A, [D1, d * D2]);
[U, S, B] = svd2(A);
if norm(diag(S)) ~= 1
S = S./norm(diag(S)); % normalisation necessary for thermal steps
end
sv = diag(S);
%% Start Truncation / Expansion
if sitej <= min(para.trustsite(end),para.L)
if sv(end) > para.svmintol
%% Expand A dims
% Expand A then do SVD since this leaves B with orthonormal
% rows. Padding after SVD would only leave zeros.
if sv(end) > para.svmaxtol
dimincratio = log10(sv(end)/para.svmaxtol)/2; % aggressive expansion if min(SV) > SVmaxtol.
elseif abs(log10(sv(end))) > 2
dimincratio = log10(sv(end)/para.svmintol)/(abs(log10(sv(end)))-1); % new scheme, also weighted by current lowest exponent, less aggressive than above
else
dimincratio = log10(sv(end)/para.svmintol)/2; % new scheme for big SV, even less aggressive than above
end
adddim = ceil(D1*dimincratio); % at least 1 for dimincratio != 0
if isfield(para,'tdvp')
if sitej ~= 1 && length(para.tdvp.maxBondDim) > 1
adddim = min(adddim, para.tdvp.maxBondDim(2) - D1);
else
adddim = min(adddim, para.tdvp.maxBondDim(1) - D1);
end
end
A = cat(1,A,zeros(adddim, d * D2));
[U, S, B] = svd2(A); % better since B contains no zeros but orthonormal vectors!
if norm(diag(S)) ~= 1
S = S./norm(diag(S)); % normalisation necessary for thermal steps
end
U = U(1:D1,:); % form old Bond Dim in columns == U=cat(2,U,addmat);
else
%% Truncate A dims
[U, S, B, err] = truncate(U,S,B);
% If Bond Dim truncated then sv(end) < para.svmaxtol
end
end
DB = size(S, 1);
if sitej ~= 1
para.D(sitej-1) = DB;
end
B = reshape(B, [DB, d, D2]); B = permute(B, [1, 3, 2]);
U = U * S;
else
error('Have not implemented different parity yet!');
end
resultsSite = sitej-1;
end
sv = diag(S);
vNE = vonNeumannEntropy(S);
if nargin == 4 && resultsSite >= 1
results.Amat_vNE(resultsSite) = vNE;
sv(sv < eps) = 0; % need to set to zero, since otherwise might have artifacts after expansion!
if length(sv)==para.D(resultsSite) || (length(sv)==para.D(resultsSite)/2 && (para.parity~='n'))
results.Amat_sv{resultsSite} = sv;
end
results.Amat_truncErr(resultsSite) = err;
end
function [U, S, V, err] = truncate(U,S,V)
%% Truncates Bond Dimensions
% keeps always one SV in range [para.svmaxtol, para.svmintol]
% if not possible -> Do nothing, expand later
% keeps always at least para.Dmin many dimensions
% by Florian Schroeder 29/10/2014
sv = diag(S);
err = 0;
%% Truncate A dims
keepdims = find(sv > para.svmintol);
% If smallest SV too large, keep 1 more
while (sv(keepdims(end))/norm(sv(keepdims)) > para.svmaxtol)
if keepdims(end)+1 <= length(sv)
keepdims = [keepdims;keepdims(end)+1];
else
% newBondDim = length(sv); % = old Dim
% expand Bond dims! Happens automatically!
end
end
if length(keepdims) < para.Dmin % keep at least Dmin bonds
keepdims = 1:para.Dmin;
end
if length(keepdims) < length(sv) % && sitej < para.trustsite(end)
U = U(:,keepdims); % keep columns
S = diag(sv(keepdims)./norm(sv(keepdims)));
V = V(keepdims,:); % keep rows
% newBondDim = length(keepdims);
% truncation error estimation
truncdims = 1:length(sv);
truncdims(keepdims) = [];
err = sum(sv(truncdims).^2); % all disregarded SV^2
end
end
end