1. 数学模型

观测数学模型可由下面公式给出
$$y=|Ax{|}^2$$
其中$x\in {\mathbb{C}}^n$， $A\in {\mathbb{C}}^{m\times n}$， $y\in {\mathbb{R}}^m$。

所以我们要求解的问题可归为如下非凸最小二乘问题
$$\underset{z\in {\mathbb{C}}^n}{\min }f(z)=\frac{1}{2}\Vert |Az{|}^2-y{\Vert }_2^2$$

2. Wirtinger Flow 算法

该算法可以总结成两步：1. 光谱初始化 2. Wirtinger梯度下降

2.1. 光谱初始化方法

具体步骤如下：

1. 能量标定系数
   $${\lambda }^2=n\frac{1_m^{\top }y}{\Vert A{\Vert }_F^2}$$

2. 构造自伴矩阵
   $$Y=\frac{1}{m}{A}^{\ast }\mathrm{diag}(y)A$$
   求得其最大特征向量$v$

3. 缩放得到初始点
   $$z_0=\lambda v$$

2.2. Wirtinger梯度下降

更新公式如下：
$$z_{\tau +1}=z_{\tau }-\frac{{\mu }_{\tau +1}}{\Vert z_0{\Vert }_2^2}\underset{\nabla f(z_{\tau })}{\underbrace{\left(\frac{1}{m}{A}^{\ast }\right[(|A{z}_{\tau }{|}^2-y)\odot (A{z}_{\tau })])}}$$
公式中的$\mu$更新根据经验公式
$${\mu }_{\tau }=\min (1-\exp (-\tau /{\tau }_0),0.2),{\tau }_0\approx 330$$

3. 算法实现

3.1. Matlab实现

clear; close all; clc
%% Measurement model 
% Signal length
n = 128; 
% Complex signal
x = randn(n,1) + 1i*randn(n,1);                     % measurement number  
m = 5 * n; 
% Measurement matrix
A = 1/sqrt(2)*randn(m,n) + 1i/sqrt(2)*randn(m,n); 
% Measured values
y = abs(A*x).^2 ;                                   %% Initialization
% power method to get the initial guess
npower_iter = 50;% Scaled coefficient lambda
lam = sqrt(n * sum(y) / norm(A, 'fro')^2);% Random input
z0 = randn(n,1); z0 = z0/norm(z0,'fro');
for tt = 1:npower_iterz0 = 1/m * A'*(y .* (A*z0));z0 = z0/norm(z0,'fro');
end% Initialized ouput
z = lam * z0;%% Gradient update
% Max number of iterations
max_iter = 2500;% update mu
tau0 = 330;                         
mu = @(t) min(1-exp(-t/tau0), 0.2); % Store relative errors
relative_error = zeros(max_iter, 1);for tt = 1:max_iter  Az = A*z;% Wirtinger gradientgrad  = 1/m* A'*( ( abs(Az).^2-y ) .* Az ); % ||z0||=lamz = z - mu(tt)/lam^2 * grad;            % Calculate relative error valuerelative_error_val = norm(x - exp(-1i*angle(trace(x'*z))) * z, 'fro')/norm(x,'fro');relative_error(tt) = relative_error_val;  
end
%%
figure,semilogy(relative_error,'LineWidth',1.8, 'Color',[0 0.4470 0.7410])
xlabel('Iteration','FontSize',16,'FontName','Times New Roman')
ylabel('Relative error','FontSize',16,'FontName','Times New Roman')
title('Wirtinger Flow Convergence','FontSize',16,'FontWeight','bold')

3.2. Python实现

import numpy as np
import matplotlib.pyplot as pltn = 128                                        # 信号长度
x = np.random.randn(n) + 1j * np.random.randn(n)   # 复值真信号m = 5 * n                                      # 测量数
A = (np.random.randn(m, n) + 1j * np.random.randn(m, n)) / np.sqrt(2)y = np.abs(A @ x) ** 2                         # 强度观测 |Ax|^2# ---------- Initialization (power method) ------------------------------------
npower_iter = 50                               # 幂迭代次数
lam = np.sqrt(n * y.sum() / np.linalg.norm(A, "fro") ** 2)  # λz0 = np.random.randn(n) + 1j * np.random.randn(n)
z0 /= np.linalg.norm(z0)for _ in range(npower_iter):z0 = (A.conj().T @ (y * (A @ z0))) / mz0 /= np.linalg.norm(z0)z = lam * z0                                   # 初值# ---------- Gradient update ---------------------------------------------------
max_iter = 2500
tau0 = 330.0
rel_err = np.zeros(max_iter)for tt in range(max_iter):mu = min(1.0 - np.exp(-(tt + 1) / tau0), 0.2)   # 步长 μ_tAz = A @ zgrad = (A.conj().T @ ((np.abs(Az) ** 2 - y) * Az)) / mz = z - (mu / lam ** 2) * grad# 相对误差theta = np.angle(np.vdot(x, z))           # vdot = x* · zrel_err[tt] = np.linalg.norm(x - np.exp(-1j * theta) * z) / np.linalg.norm(x)# ---------- plot -------------------------------------------------------
plt.figure(figsize=(6.2, 4.2), facecolor="w")
plt.semilogy(rel_err, lw=1.8, color=(0.0, 0.447, 0.741))
plt.xlabel("Iteration", fontsize=13)
plt.ylabel("Relative error", fontsize=13)
plt.title("Wirtinger Flow Convergence (1-D)", fontsize=15, weight="bold")
plt.grid(ls="--", alpha=0.3)
plt.tight_layout()
plt.show()

参考文献

Candes E J, Li X, Soltanolkotabi M. Phase retrieval via Wirtinger flow: Theory and algorithms[J]. IEEE Transactions on Information Theory, 2015, 61(4): 1985-2007.