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

1. 简介

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

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

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

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

2. 算法描述

2.1 初始化(Initialization)

  • 输入:

    • 观测数据 { y r } r = 1 m \{y_r\}_{r=1}^m {yr}r=1m, 其中 y r = ∣ ⟨ a r , x ⟩ ∣ 2 y_r=|\langle a_r,x\rangle|^2 yr=ar,x2, 采样向量 { a r } r = 1 m \{a_r\}_{r=1}^m {ar}r=1m
  • 步骤:

1.计算参数 λ : \lambda: λ:

λ 2 = n ∑ r = 1 m y r ∑ r = 1 m ∥ a r ∥ 2 \lambda^2=\frac{n\sum_{r=1}^my_r}{\sum_{r=1}^m\|a_r\|^2} λ2=r=1mar2nr=1myr

​2.构造矩阵 Y : Y: Y:
Y = 1 m ∑ r = 1 m y r a r a r ∗ Y=\frac1m\sum_{r=1}^my_ra_ra_r^* Y=m1r=1myrarar

  1. 计算 Y Y Y 的最大特征值对应的特征向量 z 0 z_0 z0, 并将其归一化为 ∥ z 0 ∥ = λ . \|z_0\|=\lambda. z0=λ.

2.2 迭代更新(Iterative Updates)

  • 输入:

    • 初始猜测 z 0 z_{0} z0
    • 学习率 μ \mu μ
    • 最大迭代次数 T T T
  • 步骤:

​ 1.设置初始值 z = z 0 z=z_0 z=z0

​ 2.对于每次迭代 τ = 0 , 1 , 2 , … , T − 1 : \tau=0,1,2,\ldots,T-1: τ=0,1,2,,T1:

​ 1.计算梯度 ∇ f ( z ) : \nabla f(z): f(z):
∇ f ( z ) = 1 m ∑ r = 1 m ( ∣ ⟨ a r , z ⟩ ∣ 2 − y r ) a r a r ∗ z \nabla f(z)=\frac{1}{m}\sum_{r=1}^{m}(|\langle a_r,z\rangle|^2-y_r)a_ra_r^*z f(z)=m1r=1m(ar,z2yr)ararz
​ 2.更新 z : z: z:
z = z − μ ∇ f ( z ) z=z-\mu\nabla f(z) z=zμf(z)

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

2.3 学习率调整(Learning Rate Adjustment)

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

3. 代码实现

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

  • 29
    点赞
  • 26
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

xy_optics

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值