Returns a manifold struct to optimize over phases of fft's of real signals
function M = realphasefactory(n)
function M = realphasefactory(n, z0)
function M = realphasefactory(n, z0, zmax)
If x is a real vector of length n, then y = fft(x) is a complex vector
which obeys certain symmetries. Specifically, for any integer k,
y(1+mod(k, n)) = conj(y(1+mod(n-k, n)))
The same holds for the phases of the Fourier transform z = sign(y).
This factory returns a Manopt manifold structure which represents the set
of complex vectors z of length n which could be the phases of the Fourier
transform of a real signal of length n:
abs(z) = 1 and z(1+mod(k, n)) = conj(z(1+mod(n-k, n))) for each k.
For k = 1, this readily implies that z(1) is +1 or -1, so that the set of
possible z's is disconnected. To choose which connected component to work
with, set the second input z0 to +1 or -1 (this is the sign of the mean
of x). By default, z0 = 1.
Furthermore, if n is even, then k = n/2 implies z(1+n/2) is +1 or -1 as
well, thus further disconnecting the set of acceptable z's. To choose
which component to work with, set the third input zmax to +1 or -1. By
default, it is +1.
The Riemannian manifold structure is the Riemannian submanifold
structure from the embedding space R^2 x ... x R^2, i.e., the complex
circles are identified with the unit circle in the real plane.
Concretely, this means the inner product is <u, v>_z = real(u'*v).
Tangent vectors at z are complex vectors of length n which notably
satisfy z(1+0) = 0 and, if n is even, z(1+n/2) = 0.
n must be integer and n >= 3 (for n = 1:2 the manifold has dimension 0).
Extra functions available in M include M.up, M.down and M.downup. They
allow to capture the symmetries concisely, as:
M.up(z) == conj(M.down(z)).
See in code for more details.
See also complexcirclefactory0001 function M = realphasefactory(n, z0, zmax) 0002 % Returns a manifold struct to optimize over phases of fft's of real signals 0003 % 0004 % function M = realphasefactory(n) 0005 % function M = realphasefactory(n, z0) 0006 % function M = realphasefactory(n, z0, zmax) 0007 % 0008 % If x is a real vector of length n, then y = fft(x) is a complex vector 0009 % which obeys certain symmetries. Specifically, for any integer k, 0010 % 0011 % y(1+mod(k, n)) = conj(y(1+mod(n-k, n))) 0012 % 0013 % The same holds for the phases of the Fourier transform z = sign(y). 0014 % 0015 % This factory returns a Manopt manifold structure which represents the set 0016 % of complex vectors z of length n which could be the phases of the Fourier 0017 % transform of a real signal of length n: 0018 % 0019 % abs(z) = 1 and z(1+mod(k, n)) = conj(z(1+mod(n-k, n))) for each k. 0020 % 0021 % For k = 1, this readily implies that z(1) is +1 or -1, so that the set of 0022 % possible z's is disconnected. To choose which connected component to work 0023 % with, set the second input z0 to +1 or -1 (this is the sign of the mean 0024 % of x). By default, z0 = 1. 0025 % 0026 % Furthermore, if n is even, then k = n/2 implies z(1+n/2) is +1 or -1 as 0027 % well, thus further disconnecting the set of acceptable z's. To choose 0028 % which component to work with, set the third input zmax to +1 or -1. By 0029 % default, it is +1. 0030 % 0031 % The Riemannian manifold structure is the Riemannian submanifold 0032 % structure from the embedding space R^2 x ... x R^2, i.e., the complex 0033 % circles are identified with the unit circle in the real plane. 0034 % Concretely, this means the inner product is <u, v>_z = real(u'*v). 0035 % Tangent vectors at z are complex vectors of length n which notably 0036 % satisfy z(1+0) = 0 and, if n is even, z(1+n/2) = 0. 0037 % 0038 % n must be integer and n >= 3 (for n = 1:2 the manifold has dimension 0). 0039 % 0040 % Extra functions available in M include M.up, M.down and M.downup. They 0041 % allow to capture the symmetries concisely, as: 0042 % 0043 % M.up(z) == conj(M.down(z)). 0044 % 0045 % See in code for more details. 0046 % 0047 % See also complexcirclefactory 0048 0049 % This file is part of Manopt: www.manopt.org. 0050 % Original author: Nicolas Boumal, Feb. 2, 2017. 0051 % Contributors: joint work with Tamir Bendory, Zhizhen Zhao and Amit Singer 0052 % Change log: 0053 % 0054 % July 20, 2017 (NB) 0055 % The distance function is now more accurate. Improved logarithm 0056 % accordingly. 0057 0058 assert(n == round(n) && n >= 3, 'n must be an integer >= 3.'); 0059 0060 even_n = (round(n/2) == n/2); 0061 0062 if ~exist('z0', 'var') || isempty(z0) 0063 z0 = 1; 0064 end 0065 if ~exist('zmax', 'var') || isempty(zmax) 0066 zmax = 1; 0067 end 0068 0069 assert(z0 == 1 || z0 == -1, 'z0 must be +1 or -1.'); 0070 assert(zmax == 1 || zmax == -1, 'zmax must be +1 or -1.'); 0071 0072 if even_n 0073 M.name = @() sprintf('Phases of fft''s of real signals of length %d (z0 = %d, zmax = %d)', n, z0, zmax); 0074 else 0075 M.name = @() sprintf('Phases of fft''s of real signals of length %d (z0 = %d)', n, z0); 0076 end 0077 0078 M.dim = @() floor((n-1)/2); 0079 0080 M.inner = @(z, v, w) real(v'*w); 0081 0082 M.norm = @(z, u) norm(u); 0083 0084 M.dist = @(z1, z2) norm(real(2*asin(.5*abs(z1 - z2)))); 0085 0086 M.typicaldist = @() pi*sqrt(n/2); 0087 0088 % Special functions to ease working with the symmetries. 0089 down = @(u) u; 0090 up = @(u) u([1 ; (n:-1:2)']); 0091 downup = @(u) (down(u) + conj(up(u)))/2; 0092 M.down = down; 0093 M.up = up; 0094 M.downup = downup; 0095 0096 M.proj = @proj; 0097 function pu = proj(z, u) 0098 duu = downup(u); 0099 pu = duu - real(duu .* conj(z)).*z; 0100 % Note that there is no need to enforce pu(1) = 0 or (if n is even) 0101 % pu(1+n/2) = 0 manually, since the IEEE standard ensures that the 0102 % above operation will be exact for those entries provided z(1) 0103 % (and possibly z(1+n/2) is +1 or -1, as should be the case. 0104 end 0105 0106 M.tangent = M.proj; 0107 0108 % For Riemannian submanifolds, converting a Euclidean gradient into a 0109 % Riemannian gradient amounts to an orthogonal projection. 0110 M.egrad2rgrad = M.proj; 0111 0112 M.ehess2rhess = @ehess2rhess; 0113 function rhess = ehess2rhess(z, egrad, ehess, zdot) 0114 rhess = M.proj(z, ehess - real(downup(egrad) .* conj(z)).*zdot); 0115 end 0116 0117 M.exp = @exponential; 0118 function y = exponential(z, v, t) 0119 0120 if nargin == 2 0121 % t = 1; 0122 tv = v; 0123 else 0124 tv = t*v; 0125 end 0126 0127 y = zeros(n, 1); 0128 0129 nrm_tv = abs(tv); 0130 0131 % We need to be careful for zero steps. 0132 mask = (nrm_tv > 0); 0133 y(mask) = z(mask).*cos(nrm_tv(mask)) + ... 0134 tv(mask).*(sin(nrm_tv(mask))./nrm_tv(mask)); 0135 y(~mask) = z(~mask); 0136 0137 end 0138 0139 M.retr = @retraction; 0140 function y = retraction(z, v, t) 0141 if nargin == 2 0142 % t = 1; 0143 tv = v; 0144 else 0145 tv = t*v; 0146 end 0147 y = sign(z+tv); 0148 end 0149 0150 M.log = @logarithm; 0151 function v = logarithm(x1, x2) 0152 v = M.proj(x1, x2 - x1); 0153 di = real(2*asin(.5*abs(x1 - x2))); 0154 nv = abs(v); 0155 factors = di ./ nv; 0156 factors(di <= 1e-6) = 1; 0157 v = v .* factors; 0158 end 0159 0160 M.hash = @(z) ['z' hashmd5( [real(z(:)) ; imag(z(:))] ) ]; 0161 0162 M.rand = @random; 0163 function z = random() 0164 z = sign(downup(randn(n, 1) + 1i*randn(n, 1))); 0165 z(1) = z0; 0166 if even_n 0167 z(1 + n/2) = zmax; 0168 end 0169 end 0170 0171 M.randvec = @randomvec; 0172 function v = randomvec(z) 0173 v = M.proj(z, randn(n, 1) + 1i*randn(n, 1)); 0174 v = v / norm(v); 0175 end 0176 0177 M.lincomb = @matrixlincomb; 0178 0179 M.zerovec = @(z) zeros(n, 1); 0180 0181 M.transp = @(z1, z2, u) M.proj(z2, u); 0182 0183 M.pairmean = @pairmean; 0184 function z = pairmean(z1, z2) 0185 z = sign(z1+z2); 0186 end 0187 0188 % This vec/mat pair is an isometry which allows to switch between the 0189 % classical representation of tangent vectors---as complex vectors of 0190 % length n---to real vectors of length M.dim() whose entries are the 0191 % coordinates of the tangent vector in the basis 1i*z, for the first 0192 % half. A scaling of sqrt(2) is applied to ensure isometry, since 0193 % tangent vectors are represented with only half of their entries. 0194 I = 2 : floor((n+1)/2); 0195 if even_n 0196 middle = 0; 0197 else 0198 middle = []; 0199 end 0200 M.vec = @(z, u_mat) sqrt(2)*real(u_mat(I) .* conj(1i*z(I))); 0201 M.mat = @(z, u_vec) [0 ; u_vec.*(1i*z(I)) ; middle ; ... 0202 flipud(conj(u_vec.*(1i*z(I))))]/sqrt(2); 0203 M.vecmatareisometries = @() true; 0204 0205 end