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

文章目录

    • 1. 简介
    • 2. 算法描述
      • 2.1 初始化(Initialization)
      • 2.2 迭代更新(Iterative Updates)
      • 2.3 学习率调整(Learning Rate Adjustment)
    • 3. 代码实现
      • 3.1 一维信号测试 (Gaussian model)
      • 3.2 一维信号测试 (Coded diffraction model)
      • 3.2 二维图片测试
    • 参考文献

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 @ zgrad = (1 / m) * A.conj().T @ ((np.abs(yz)**2 - y) * yz)z = z - mu(t) / normest**2 * gradRelerrs.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) * Bzgrad = At(C)z = z - mu(t) / normest**2 * gradRelerrs.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) * Bzgrad = At(C)z = z - mu(t) / normest**2 * gradRelerrs.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.

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.rhkb.cn/news/335964.html

如若内容造成侵权/违法违规/事实不符,请联系长河编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

最长递增子序列,交错字符串

第一题&#xff1a; 代码如下&#xff1a; int lengthOfLIS(vector<int>& nums) {//dp[i]表示以第i个元素为结尾的最长子序列的长度int n nums.size();int res 1;vector<int> dp(n, 1);for (int i 1; i < n; i){for (int j 0; j < i; j){if (nums[i]…

ArcGIS提取含有计曲线的等高线

喜欢就关注我们吧&#xff01; 今天我么来看看&#xff0c;如何利用DEM提取含有计曲线的等高线&#xff01; 常规的话我们利用DEM提取的等高线都是不带计曲线的&#xff0c;无法把计曲线标注出来&#xff0c;今天我们就来看下&#xff0c;如何处理一下哦&#xff01;提取带有计…

Python内置方法串讲:类型转化与实用技巧

新书上架~&#x1f447;全国包邮奥~ python实用小工具开发教程http://pythontoolsteach.com/3 欢迎关注我&#x1f446;&#xff0c;收藏下次不迷路┗|&#xff40;O′|┛ 嗷~~ 目录 一、类型转化&#xff1a;从A到B的魔法 二、实用技巧&#xff1a;避免类型错误 三、总结 一…

【深度好文】AI企业融合联盟营销,做的好就是最大赢家!

AI工具市场正在迅速发展&#xff0c;现仍有不少企业陆续涌出&#xff0c;那么如何让你的工具受到目标群体的关注呢&#xff1f;这相比是AI工具营销人员一直在思考的问题。 即使这个市场正蓬勃发展&#xff0c;也无法保证营销就能轻易成功。AI工具虽然被越来越多人认可和接受&a…

JAVA类与方法·易错题分析

分析一下作业中关于类与方法写错或者易错的题。 N o . 1 No.1 No.1 下面程序的执行结果是______。 public class Test7 {public static void main(String[] args){new B().display();} } class A{public void draw() {System.out.print("Draw A.");}public void di…

金蝶云星空数据库迁移后,显示 error: 40 - 无法打开到 SQL Server 的连接的解决方法

原因&#xff1a;数据库迁移/或者更新IP后&#xff0c;与之前添加的数据库地址不一致导致无法连接数据库&#xff1b; 解决方法&#xff1a;修改IP为目前数据库的IP&#xff1b; 文件路径&#xff1a;在ManageSite\APP_Data\Common.config中&#xff0c;修改DbServerInstance…

动态规划part03 Day43

LC343整数拆分&#xff08;未掌握&#xff09; 未掌握分析&#xff1a;dp数组的含义没有想清楚&#xff0c;dp[i]表示分解i能够达到的最大乘积&#xff0c;i能够如何分解呢&#xff0c;从1开始遍历&#xff0c;直到i-1&#xff1b;每次要不是j和i-j两个数&#xff0c;要不是j和…

Function Calling学习

Function Calling第一篇 Agent&#xff1a;AI 主动提要求Function Calling&#xff1a;AI 要求执行某个函数场景举例&#xff1a;明天上班是否要带伞&#xff1f;AI反过来问你&#xff0c;明天天气怎么样&#xff1f; Function Calling 的基本流程 Function Calling 完整的官…

adb的常见操作和命令

最近学习adb的时候&#xff0c;整理了一些adb的使用场景&#xff0c;如&#xff1a;adb与设备交互&#xff0c;adb的安装、卸载&#xff0c;adb命令启动&#xff0c;通过命令清除缓存&#xff0c;文件传输和日志操作。 adb的两大作用&#xff1a;在app测试的时候可以提供监控日…

CentOS8搭载正反向解析dns服务器

1.介绍&#xff08;是什么&#xff09; DNS&#xff08;Domain Name System&#xff09;&#xff0c;即域名系统&#xff0c;是一个将域名和 IP 地址相互映射的分布式数据库&#xff0c;它可以将用户输入的域名转换成对应的 IP 地址。DNS 由多个服务器组成&#xff0c;分为多个…

git 学习随笔

git 学习随笔 基本概念 git 对待数据类似快照流的形式而不是类似 cvs 那样的纪录文件随时间逐步积累的差异 git 中所有数据在存储钱都会计算校验和&#xff08;hash) 三种状态&#xff1a;已提交(committed)&#xff0c;已修改(modified)&#xff0c;已暂存(staged)。 add…

【spring】@PathVariable注解学习

PathVariable介绍 PathVariable是Spring框架中的一个注解&#xff0c;主要用于处理RESTful风格URL中的路径变量。在RESTful接口设计中&#xff0c;我们经常将资源的ID或者其他标识信息直接放在URL路径中&#xff0c;而不是作为查询参数。PathVariable注解使得控制器方法能够轻…

【Linux】Linux基本指令2

我们接着上一篇&#xff1a;http://t.csdnimg.cn/bSJx8 我们接着完善ls指令 我们可以直接匹配对应格式的文件匹配出来 1.man指令&#xff08;重要&#xff09;&#xff1a; Linux的命令有很多参数&#xff0c;我们不可能全记住&#xff0c;我们可以通过查看联机手册获取帮助…

原神抽卡点名程序教程(直接下载用)

今天我要给大家分享一个在抖音上特别火的视频——原神抽卡点名程序教程。 &#xff08;要源码的私信扣31&#xff09; 废话不多说&#xff0c;直接上效果图 &#xff1a; 步骤1&#xff1a; 步骤2&#xff1a;&#xff08;写名单&#xff0c;前面加数字代表星级&#xff0c;用…

去除字符串中的空格和特殊字符

自学python如何成为大佬(目录):https://blog.csdn.net/weixin_67859959/article/details/139049996?spm1001.2014.3001.5501 用户在输入数据时&#xff0c;可能会无意中输入多余的空格&#xff0c;或在一些情况下&#xff0c;字符串前后不允许出现空格和特殊字符&#xff0c;…

【C语言】自定义类型:联合体和枚举

1. 联合体 1.1 联合体的特点 像结构体⼀样&#xff0c;联合体也是由⼀个或者多个成员构成&#xff0c;这些成员可以不同的类型。 但是编译器只为最⼤的成员分配⾜够的内存空间。联合体的特点是所有成员共⽤同⼀块内存空间所以联合体也叫&#xff1a;共⽤体。 union Un {char…

DETR整体模型结构解析

DETR流程 Backbone用卷积神经网络抽特征。最后通过一层1*1卷积转化到d_model维度fm&#xff08;B,d_model,HW&#xff09;。 position embedding建立跟fm维度相同的位置编码(B&#xff0c;d_model,HW&#xff09;。 Transformer Encoder,V为fm&#xff0c;K&#xff0c;Q为fm…

【Python】解决Python报错:TypeError: can only concatenate str (not “int“) to str

&#x1f9d1; 博主简介&#xff1a;阿里巴巴嵌入式技术专家&#xff0c;深耕嵌入式人工智能领域&#xff0c;具备多年的嵌入式硬件产品研发管理经验。 &#x1f4d2; 博客介绍&#xff1a;分享嵌入式开发领域的相关知识、经验、思考和感悟&#xff0c;欢迎关注。提供嵌入式方向…

反转!Greenplum 还在,快去 Fork 源码

↑ 关注“少安事务所”公众号&#xff0c;欢迎⭐收藏&#xff0c;不错过精彩内容~ 今早被一条消息刷爆群聊&#xff0c;看到知名开源数仓 Greenplum 的源码仓“删库跑路”了。 要知道 GP 新东家 Broadcom 前几日才刚刚免费开放了 VMware Workstation PRO 17 和 VMware Fusion P…

Selenium 自动化测试工具(1) (Selenium 工作原理,常用API的使用)

文章目录 什么是自动化测试什么是测试工具&#xff1a;Selenium 工作原理(重要)Selenium API定位元素CSS 选择器xpath 定位元素 通过Java代码实现自动化1. 定位元素2. 关闭浏览器3. 获取元素文本4. 鼠标点击与键盘输入5. 清空内容6.打印信息 什么是自动化测试 关于自动化&…