通过Wirtinger流进行相位恢复：理论与算法

1. 简介

Wirtinger流，包含通过谱方法初始化并使用类似梯度下降的新更新规则迭代地细化估计。这种方法承诺通过最少数量的随机测量实现精确的相位恢复，同时在计算和数据资源方面都具有高效性。

传统方法，如Gerchberg-Saxton算法及其变体，通过迭代地应用投影，使频谱的幅度与观察数据匹配。这些方法依赖于关于信号的先验知识，并且缺乏理论上的收敛保证。

Wirtinger流算法包括两个主要步骤：

1. 通过谱方法初始化: 初始猜测 $z_o$ 是通过计算从测量数据构建的矩阵的主特征向量得到的。
2. 类似梯度下降的更新: 使用非凸目标函数的梯度迭代地细化 $z$。步长动态调整以确保收敛。

2. 算法描述

2.1 初始化(Initialization)

• 输入：

  • 观测数据 { y r } r = 1 m \{y_r\}_{r=1}^m {yr​}r=1m​, 其中 $y_r=|⟨a_r,x⟩{|}^2$， 采样向量 { a r } r = 1 m \{a_r\}_{r=1}^m {ar​}r=1m​

• 步骤：

1.计算参数 $\lambda :$​

$${\lambda }^2=\frac{n{\sum }_{r=1}^m{y}_r}{{\sum }_{r=1}^m\Vert a_r{\Vert }^2}$$

​2.构造矩阵 $Y:$
$$Y=\frac{1}{m}\sum _{r=1}^m{y}_r{a}_r{a}_r^{\ast }$$

3. 计算 $Y$ 的最大特征值对应的特征向量 $z_0$, 并将其归一化为 $\Vert z_0\Vert =\lambda .$

2.2 迭代更新(Iterative Updates)

• 输入：

  • 初始猜测 $z_0$
  • 学习率 $\mu$
  • 最大迭代次数 $T$

• 步骤：

​ 1.设置初始值 $z=z_0$

​ 2.对于每次迭代 $\tau =0,1,2,\dots ,T-1:$

​ 1.计算梯度 $\nabla f(z):$
$$\nabla f(z)=\frac{1}{m}\sum _{r=1}^m(|⟨a_r,z⟩{|}^2-y_r)a_r{a}_r^{\ast }z$$
​ 2.更新 $z:$
$$z=z-\mu \nabla f(z)$$

• 输出：
  • 恢复的信号 $z$

2.3 学习率调整（Learning Rate Adjustment）

​ 在迭代过程中，可以根据梯度估计的不确定性调整学习率 $\mu$。通常，早期迭代使用较小的学习率，随后逐渐增大。建议的学习率更新公式为：
$${\mu }_{\tau }=\min (1-e^{-\tau /{\tau }_0},{\mu }_{\max })$$
其中${\tau }_0$ 和 ${\mu }_{max}$ 是实验中调整的超参数。

3. 代码实现

1. 文章中的方法：
   • 通过构造矩阵 $Y$ 并计算其最大特征值对应的特征向量，来获得初始估计。
   • 这种方法直接利用特征值分解，得到了与最大特征值对应的特征向量。
2. 程序中的方法：
   • 使用幂法迭代，这是一种求矩阵最大特征值和特征向量的数值方法。
   • 幂法迭代通过反复应用线性算子 $A$和其共轭转置 $At$，逐步逼近最大特征值对应的特征向量。

3.1 一维信号测试 （Gaussian model）

import numpy as np
import matplotlib.pyplot as plt

# 生成信号和数据
n = 128
x = np.random.randn(n) + 1j * np.random.randn(n)

m = round(4.5 * n)
A = (1 / np.sqrt(2)) * (np.random.randn(m, n) + 1j * np.random.randn(m, n))
y = np.abs(A @ x)**2

# 初始化
npower_iter = 50
z0 = np.random.randn(n) + 1j * np.random.randn(n)
z0 = z0 / np.linalg.norm(z0)

for _ in range(npower_iter):
    z0 = A.conj().T @ (y * (A @ z0))
    z0 = z0 / np.linalg.norm(z0)

normest = np.sqrt(np.sum(y) / len(y))
z = normest * z0
Relerrs = [np.linalg.norm(x - np.exp(-1j * np.angle(np.vdot(x, z))) * z) / np.linalg.norm(x)]

# 主循环
T = 2500
tau0 = 330
mu = lambda t: min(1 - np.exp(-t / tau0), 0.2)

for t in range(1, T + 1):
    yz = A @ z
    grad = (1 / m) * A.conj().T @ ((np.abs(yz)**2 - y) * yz)
    z = z - mu(t) / normest**2 * grad
    Relerrs.append(np.linalg.norm(x - np.exp(-1j * np.angle(np.vdot(x, z))) * z) / np.linalg.norm(x))

# 检查结果
print(f'Relative error after initialization: {Relerrs[0]:.6f}')
print(f'Relative error after {T} iterations: {Relerrs[-1]:.6f}')

# 绘制相对误差图
plt.figure()
plt.semilogy(range(T + 1), Relerrs)
plt.xlabel('Iteration')
plt.ylabel('Relative error (log10)')
plt.title('Relative error vs. iteration count')
plt.show()

3.2 一维信号测试 （Coded diffraction model）

import numpy as np
import matplotlib.pyplot as plt

# 生成信号
n = 128
x = np.random.randn(n) + 1j * np.random.randn(n)

# 生成掩码和线性采样算子
L = 6  # 掩码数量

# 生成随机相位掩码
Masks = np.random.choice([1j, -1j, 1, -1], size=(n, L))

# 随机缩放
temp = np.random.rand(n, L)
Masks = Masks * ((temp <= 0.2) * np.sqrt(3) + (temp > 0.2) / np.sqrt(2))

# 定义线性采样算子
def A(I):
    return np.fft.fft(np.conj(Masks) * np.tile(I[:, np.newaxis], (1, L)), axis=0)

def At(Y):
    return np.mean(Masks * np.fft.ifft(Y, axis=0), axis=1)

# 测量数据
Y = np.abs(A(x))**2

# 初始化
npower_iter = 50
z0 = np.random.randn(n) + 1j * np.random.randn(n)
z0 = z0 / np.linalg.norm(z0)

for _ in range(npower_iter):
    z0 = At(Y * A(z0))
    z0 = z0 / np.linalg.norm(z0)

normest = np.sqrt(np.sum(Y) / Y.size)
z = normest * z0
Relerrs = [np.linalg.norm(x - np.exp(-1j * np.angle(np.vdot(x, z))) * z) / np.linalg.norm(x)]

# 主循环
T = 2500
tau0 = 330
mu = lambda t: min(1 - np.exp(-t / tau0), 0.2)

for t in range(1, T + 1):
    Bz = A(z)
    C = (np.abs(Bz)**2 - Y) * Bz
    grad = At(C)
    z = z - mu(t) / normest**2 * grad
    Relerrs.append(np.linalg.norm(x - np.exp(-1j * np.angle(np.vdot(x, z))) * z) / np.linalg.norm(x))

# 检查结果
print(f'Relative error after initialization: {Relerrs[0]:.6f}')
print(f'Relative error after {T} iterations: {Relerrs[-1]:.6f}')

# 绘制相对误差图
plt.figure()
plt.semilogy(range(T + 1), Relerrs)
plt.xlabel('Iteration')
plt.ylabel('Relative error (log10)')
plt.title('Relative error vs. iteration count')
plt.show()

3.2 二维图片测试

import numpy as np
import matplotlib.pyplot as plt
import time
import cv2

# 加载图像
# 使用 OpenCV 读取本地图像
amplitude_image = cv2.imread("./standard_test_images/cameraman.tif", cv2.IMREAD_GRAYSCALE)
phase_image = cv2.imread("./standard_test_images/lena_gray_256.tif", cv2.IMREAD_GRAYSCALE)

# 调整图像大小
img_size = 256
amplitude_image = cv2.resize(amplitude_image, (img_size, img_size))
phase_image = cv2.resize(phase_image, (img_size, img_size))

# 将图像归一化到 [0, 1] 范围内
amplitude_image = amplitude_image / np.max(amplitude_image)
phase_image = phase_image / np.max(phase_image)

# 将振幅图像和相位图像结合生成复数图像 x
x = amplitude_image * np.exp(1j * 2 * np.pi * phase_image)

# 生成掩码和线性采样算子
L = 10  # 掩码数量
n1, n2 = amplitude_image.shape
Masks = np.zeros((n1, n2, L), dtype=complex)

# 生成随机相位掩码
for ll in range(L):
    Masks[:, :, ll] = np.random.choice([1j, -1j, 1, -1], size=(n1, n2))

# 随机缩放
temp = np.random.rand(n1, n2, L)
Masks *= (temp <= 0.2) * np.sqrt(3) + (temp > 0.2) / np.sqrt(2)

# 定义线性采样算子
def A(I):
    return np.fft.fft2(np.conj(Masks) * np.tile(I[:, :, np.newaxis], (1, 1, L)), axes=(0, 1))

def At(Y):
    return np.mean(Masks * np.fft.ifft2(Y, axes=(0, 1)), axis=2)

# 测量数据
Y = np.abs(A(x))**2


# 初始化
npower_iter = 50
z0 = np.random.randn(n1, n2)
z0 = z0 / np.linalg.norm(z0, 'fro')

# 幂法迭代
start_time = time.time()
for _ in range(npower_iter):
    z0 = At(Y * A(z0))
    z0 = z0 / np.linalg.norm(z0, 'fro')
init_time = time.time() - start_time

# 估计标准差并缩放特征向量
normest = np.sqrt(np.sum(Y) / Y.size)
z = normest * z0

# 初始相对误差
Relerrs = [np.linalg.norm(x - np.exp(-1j * np.angle(np.trace(np.conj(x).T @ z))) * z, 'fro') / np.linalg.norm(x, 'fro')]
Times = [init_time]

# 迭代优化
T = 500
tau0 = 330
mu = lambda t: min(1 - np.exp(-t / tau0), 0.4)

start_time = time.time()
for t in range(T):
    Bz = A(z)
    C = (np.abs(Bz)**2 - Y) * Bz
    grad = At(C)
    z = z - mu(t) / normest**2 * grad
    Relerrs.append(np.linalg.norm(x - np.exp(-1j * np.angle(np.trace(np.conj(x).T @ z))) * z, 'fro') / np.linalg.norm(x, 'fro'))
    Times.append(time.time() - start_time)

rec_amplitude = np.abs(z)
rec_phase = np.angle(z)

# 输出结果
print('All done!')

# 结果检查
print(f'Run time of initialization: {Times[0]:.2f}')
print(f'Relative error after initialization: {Relerrs[0]:.6f}')
print('\n')
print(f'Loop run time: {Times[-1]:.2f}')
print(f'Relative error after {T} iterations: {Relerrs[-1]:.6f}')

# 绘制相对误差图
plt.figure()
plt.semilogy(range(T + 1), Relerrs)
plt.xlabel('Iteration')
plt.ylabel('Relative error (log10)')
plt.title('Relative error vs. iteration count')
plt.show()

# 设置绘图的字体和字号
font = {'family': 'serif',
        'weight': 'normal',
        'size': 20}
plt.rc('font', **font)

# 显示结果
plt.figure(figsize=(12, 6))

plt.subplot(2, 2, 1)
plt.title('Original Amplitude')
plt.imshow(amplitude_image, cmap='gray')
plt.axis('off')

plt.subplot(2, 2, 2)
plt.title('Original Phase')
plt.imshow(phase_image, cmap='gray')
plt.axis('off')

plt.subplot(2, 2, 3)
plt.title('Rec Amplitude')
plt.imshow(rec_amplitude, cmap='gray')
plt.axis('off')

plt.subplot(2, 2, 4)
plt.title('Rec Phase')
plt.imshow(rec_phase, cmap='gray')
plt.axis('off')

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.