- 张量相关总变差(t-CTV)
- 定义和性质
- t-CTV 在张量恢复中的应用
- 张量恢复模型与理论
- 张量补全(TC)模型
- 张量鲁棒主成分分析(TRPCA)模型
- 理论保证和恢复性分析
4. 张量相关总变差(t-CTV)
在张量恢复问题中,通常需要利用张量的某些先验知识来指导恢复过程。常用的先验知识包括低秩性(low-rankness)和平滑性(smoothness)。低秩性表示张量可以在某个低维子空间中较好地表示,而平滑性表示张量沿某些方向上的变化是连续且平滑的。
t-CTV 正则化项
t-CTV(Tensor Correlated Total Variation,张量相关全变差)是一种新的正则化项,它能够同时编码张量的低秩和平滑先验。t-CTV的定义如下:
∥ T ∥ t-CTV = 1 γ ∑ k ∈ Γ ∥ G k ∥ ∗ , L \|T\|_{\text{t-CTV}} = \frac{1}{\gamma} \sum_{k \in \Gamma} \|G_k\|_{*,L} ∥T∥t-CTV=γ1k∈Γ∑∥Gk∥∗,L
其中:
- G k G_k Gk 是张量 T \mathcal{T} T 沿第 k k k 模的梯度张量。
- Γ \Gamma Γ 是包含张量沿某些方向上具有平滑连续性的先验集合。
- ∥ ⋅ ∥ ∗ , L \| \cdot \|_{*,L} ∥⋅∥∗,L 是在高阶 t-SVD 框架下的张量核范数。
- γ \gamma γ 是一个归一化常数。
梯度张量 G k G_k Gk
张量 T \mathcal{T} T 沿第 k k k 模的梯度张量 G k G_k Gk 定义为: G k = ∇ k ( T ) = T × k D n k G_k = \nabla_k(T) = T \times_k D_{n_k} Gk=∇k(T)=T×kDnk
其中,
- ∇ k ( T ) \nabla_k(T) ∇k(T)表示沿第 k k k 模上的梯度运算
- × k \times_k ×k 表示张量在第 k k k 模上的乘积运算
- D n k D_{n_k} Dnk 是一个行循环矩阵,表示沿第 k k k 模的差分操作。
- 梯度张量 G k G_k Gk 捕捉了张量在第 k k k 模上的变化。如果张量在第 K K K 模上变化平滑,那么 G k G_k Gk 的值将较小。
- 通过对 G k G_k Gk 施加核范数约束,t-CTV 正则化项鼓励张量在这些方向上的变化是平滑的。这是因为核范数的最小化倾向于减少梯度的非零元素,从而使得张量在这些方向上的变化更加平滑。
核范数 ∥ ⋅ ∥ ∗ , L \| \cdot \|_{*,L} ∥⋅∥∗,L
在高阶 t-SVD 框架下的张量核范数 ∥ ⋅ ∥ ∗ , L \| \cdot \|_{*,L} ∥⋅∥∗,L 定义为对张量进行 t-SVD 分解后,对得到的奇异值张量的核范数进行求和。
解释
- 低秩先验:t-CTV 中的核范数部分( ∥ G k ∥ ∗ , L \|G_k\|_{*,L} ∥Gk∥∗,L )诱导了张量的低秩性。这是因为核范数最小化倾向于使得张量的奇异值尽可能多地为零,从而达到低秩的效果。
- 平滑先验:t-CTV 中的梯度张量部分( G k G_k Gk )诱导了张量的平滑性。梯度张量捕捉了张量在某个方向上的变化,通过对其进行约束,保证了张量在这些方向上的变化是平滑的。
总结
t-CTV 正则化项通过结合低秩先验和平滑先验,提高了张量恢复问题的性能。具体来说,t-CTV能够在恢复过程中同时保持张量的低秩性和沿某些方向的平滑性,从而更准确地恢复原始张量。
5. 张量恢复模型与理论
本文提出了基于 t-CTV 的张量补全(Tensor Completion, TC)和张量鲁棒主成分分析(Tensor Robust Principal Component Analysis, TRPCA)模型。
- 张量补全(TC):目标是通过利用张量的低秩性和平滑性从部分观测的数据中恢复完整的张量。
- 张量鲁棒主成分分析(TRPCA):目标是从观测数据中同时分离出低秩部分和稀疏噪声。
模型定义
-
张量补全(TC)模型
目标是最小化 t-CTV 正则化项,要求恢复后的张量在观测位置与原始张量一致。数学表达如下:
min T ∥ T ∥ t-CTV s.t. P Ω ( T ) = P Ω ( T 0 ) \min_T \|T\|_{\text{t-CTV}} \quad \text{s.t.} \quad P_\Omega(T) = P_\Omega(T_0) Tmin∥T∥t-CTVs.t.PΩ(T)=PΩ(T0)
其中:
- ∥ T ∥ t-CTV \|T\|_{\text{t-CTV}} ∥T∥t-CTV 是 t-CTV 正则化项,用于编码张量的低秩性和平滑性。
- P Ω P_\Omega PΩ 是观测投影算子,只作用于已观测的位置。
- T 0 T_0 T0 是原始张量(带有缺失值)。
-
张量鲁棒主成分分析(TRPCA)模型
目标是最小化 t-CTV 正则化项和稀疏噪声的 ℓ 1 \ell_1 ℓ1 范数,同时保证恢复后的低秩张量和噪声之和等于观测数据。
min T , E ∥ T ∥ t-CTV + λ ∥ E ∥ 1 s.t. M = T + E \min_{T,E} \|T\|_{\text{t-CTV}} + \lambda \|E\|_1 \quad \text{s.t.} \quad M = T + E T,Emin∥T∥t-CTV+λ∥E∥1s.t.M=T+E
其中:
- ∥ T ∥ t-CTV \|T\|_{\text{t-CTV}} ∥T∥t-CTV 是 t-CTV 正则化项,用于编码低秩性和平滑性。
- ∥ E ∥ 1 \|E\|_1 ∥E∥1 是噪声张量 E E E 的 ℓ 1 \ell_1 ℓ1 范数,用于编码稀疏性。
- λ \lambda λ 是平衡参数,用于调节低秩部分和稀疏噪声之间的权重。
- M M M 是观测数据张量。
模型的作用
- 张量补全(TC)模型:通过最小化 t-CTV 正则化项,保证了恢复后的张量既具有低秩性又在某些方向上是平滑的,从而能够更准确地填补缺失值。
- 张量鲁棒主成分分析(TRPCA)模型:通过同时最小化 t-CTV 正则化项和噪声的 ℓ 1 \ell_1 ℓ1 范数,能够有效分离出低秩结构和稀疏噪声,从而实现对观测数据的鲁棒分解。
总变差(Total Variation, TV)
总变差(Total Variation, TV)是图像处理和信号处理中常用的一种正则化技术,用于捕捉图像或信号的平滑性和边缘特性。它在图像去噪、图像修复、压缩感知等领域有着广泛的应用。
总变差的基本思想是最小化图像或信号的梯度总和,从而达到去除噪声的效果,同时保留图像的边缘和重要细节。总变差可以分为两种形式:各向异性 TV(TV-1)和各向同性 TV(TV-2)。
各向异性 TV(TV-1)
在各向异性 TV 中,总变差定义为所有像素在水平方向和垂直方向梯度的绝对值之和。具体来说,给定一个图像 I I I ,其 TV-1 可以表示为:
TV-1 ( I ) = ∑ i , j ( ∣ I i + 1 , j − I i , j ∣ + ∣ I i , j + 1 − I i , j ∣ ) \text{TV-1}(I) = \sum_{i,j} \left( |I_{i+1,j} - I_{i,j}| + |I_{i,j+1} - I_{i,j}| \right) TV-1(I)=i,j∑(∣Ii+1,j−Ii,j∣+∣Ii,j+1−Ii,j∣)
其中 I i , j I_{i,j} Ii,j 表示图像在位置 ( i , j ) (i, j) (i,j) 的像素值。
各向同性 TV(TV-2)
在各向同性 TV 中,总变差定义为所有像素在水平方向和垂直方向梯度的平方和的平方根的总和。具体来说,给定一个图像 I I I ,其 TV-2 可以表示为:
TV-2 ( I ) = ∑ i , j ( I i + 1 , j − I i , j ) 2 + ( I i , j + 1 − I i , j ) 2 \text{TV-2}(I) = \sum_{i,j} \sqrt{(I_{i+1,j} - I_{i,j})^2 + (I_{i,j+1} - I_{i,j})^2} TV-2(I)=i,j∑(Ii+1,j−Ii,j)2+(Ii,j+1−Ii,j)2
-
为什么边缘部分得到了很好的保留
总变差(TV)正则化能够很好地保留图像的边缘,是因为**它在优化过程中倾向于保留图像中的显著变化,而不是平滑所有的细节。**这里是详细的解释:
总变差的工作原理
总变差正则化的核心思想是最小化图像梯度的总和,从而使图像在大部分区域保持平滑,但在边缘区域保持显著的变化。这种特性源于以下几点:
- 梯度的定义:图像的梯度表示像素值的变化率。在平滑区域,梯度值很小,而在边缘区域,梯度值很大。通过最小化总变差,算法会倾向于减少平滑区域的梯度变化,从而去除噪声。
- 各向异性 TV 和各向同性 TV:
- 各向异性 TV (TV-1) 通过对水平方向和垂直方向的梯度分别进行计算,这样在处理噪声时可以有效减少噪声,但在边缘区域仍能保持大的梯度变化。
- 各向同性 TV (TV-2) 通过计算梯度的欧几里得范数来考虑所有方向上的变化,这样可以更自然地保留图像的几何特性和边缘。
- 边缘的保护机制:在总变差正则化中,算法会更倾向于在大部分区域使梯度变小,但在边缘区域,由于梯度值大,总变差的贡献显著,因此优化过程不会过度平滑这些区域,从而保留了边缘。
例子说明
假设我们有一个简单的图像,它有一个显著的边缘:
[100, 100, 100, 100, 0, 0, 0, 0] [100, 100, 100, 100, 0, 0, 0, 0] [100, 100, 100, 100, 0, 0, 0, 0] [100, 100, 100, 100, 0, 0, 0, 0]
在这个图像中,左侧和右侧之间有一个明显的边缘。如果我们使用总变差去噪,它会试图在平滑区域(例如左侧和右侧的均匀区域)减少梯度变化,但在边缘处(100 与 0 之间的过渡)保持显著的变化。
差分矩阵
- 前向差分矩阵:
- 后向差分矩阵:
- 中心差分矩阵:
- 周期性差分矩阵:
不同类型的差分矩阵适用于不同的应用场景,具体如下:
1. 前向差分矩阵 (Forward Difference Matrix)
D forward = [ − 1 1 0 0 0 − 1 1 0 0 0 − 1 1 0 0 0 − 1 ] D_{\text{forward}} = \begin{bmatrix} -1 & 1 & 0 & 0 \\ 0 & -1 & 1 & 0 \\ 0 & 0 & -1 & 1 \\ 0 & 0 & 0 & -1 \end{bmatrix} Dforward= −10001−10001−10001−1
应用场景:
- 时间序列分析:用于计算时间序列数据中的前向变化率,例如股票价格、传感器数据等。
- 边界条件:适用于具有明确起始边界的情况,例如计算某一段时间内的累积变化。
优点:
- 计算简单,适用于在线计算和实时数据处理。
- 对起始边界处理较好。
缺点:
- 对数据的末尾点没有充分利用,可能导致较大的误差。
2. 后向差分矩阵 (Backward Difference Matrix)
D backward = [ 1 0 0 0 − 1 1 0 0 0 − 1 1 0 0 0 − 1 1 ] D_{\text{backward}} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ -1 & 1 & 0 & 0 \\ 0 & -1 & 1 & 0 \\ 0 & 0 & -1 & 1 \end{bmatrix} Dbackward= 1−10001−10001−10001
应用场景:
- 时间序列回溯:用于计算过去时间点的变化率,例如分析历史数据的变化趋势。
- 边界条件:适用于具有明确结束边界的情况。
优点:
- 计算过去数据点的变化率,适用于回溯分析。
- 对数据的起始点没有利用,适用于末尾边界处理。
缺点:
- 对数据的起始点没有充分利用,可能导致较大的误差。
3. 中心差分矩阵 (Central Difference Matrix)
D central = 1 2 [ 0 1 0 − 1 − 1 0 1 0 0 − 1 0 1 1 0 − 1 0 ] D_{\text{central}} = \frac{1}{2} \begin{bmatrix} 0 & 1 & 0 & -1 \\ -1 & 0 & 1 & 0 \\ 0 & -1 & 0 & 1 \\ 1 & 0 & -1 & 0 \end{bmatrix} Dcentral=21 0−10110−10010−1−1010
应用场景:
- 图像处理:用于计算图像梯度,例如边缘检测、纹理分析等。
- 数值微分:在数值计算中用于逼近导数值,具有较高的精度。
优点:
- 对称性好,误差较小。
- 适用于数据均匀分布的情况,能够提供更精确的导数近似。
缺点:
- 计算复杂度稍高于前向和后向差分。
- 对边界处理较弱,需要特殊处理边界点。
4. 周期性差分矩阵 (Periodic Difference Matrix)
D periodic = [ − 1 1 0 0 0 − 1 1 0 0 0 − 1 1 1 0 0 − 1 ] D_{\text{periodic}} = \begin{bmatrix} -1 & 1 & 0 & 0 \\ 0 & -1 & 1 & 0 \\ 0 & 0 & -1 & 1 \\ 1 & 0 & 0 & -1 \end{bmatrix} Dperiodic= −10011−10001−10001−1
应用场景:
- 周期性数据:适用于处理周期性数据,例如气象数据、信号处理中的频率分析等。
- 图像处理:用于处理具有周期性边界条件的图像数据,例如拼接图像时的边缘处理。
优点:
- 能够处理周期性边界条件,适用于循环数据。
- 保证周期性数据的连续性,避免边界突变。
缺点:
- 对非周期性数据不适用,可能导致误差。
porder_diff
- 计算张量
X
在指定方向index
上的差分。 - 返回一个与
X
尺寸相同的张量,表示在该方向上的差分结果。
function DX = porder_diff(X, direction) % 计算沿指定方向的差分张量(梯度图)
% 获取张量X的尺寸
dim = size(X);
% 初始化索引,适用于多维数组
index_first = repmat({':'}, 1, ndims(X)); % 创建一个元胞数组,长度为 X 的维度数,每个元素都是':'
index_first(direction) = {1}; % 将指定方向的索引设置为1(表示第一个元素)
index_end = repmat({':'}, 1, ndims(X)); % 创建另一个元胞数组,长度为X的维度数,每个元素都是':'
index_end(direction) = {dim(direction)}; % 将指定方向的索引设置为该维度的最后一个元素
% 计算起始切片和结束切片的差
slice = X(index_first{:}) - X(index_end{:});
% 计算沿指定方向的差分
DX = diff(X, 1, direction);
% 将计算的差分和切片差拼接起来
DX = cat(direction, DX, slice);
-
详细解释:
-
获取张量
X
的尺寸dim = size(X);
dim
是一个包含X
各个维度尺寸的向量。例如,如果X
是一个 3D 张量,其尺寸为[4, 5, 6]
,那么dim
将是[4, 5, 6]
。
-
初始化索引
index_first = repmat({':'}, 1, ndims(X)); index_end = repmat({':'}, 1, ndims(X));
index_first
和index_end
初始化为适用于多维数组的全选索引。例如,对于一个 3D 张量,这些索引初始为{':', ':', ':'}
。表示全选该维度上的所有元素。
-
设置起始和结束索引
index_first(direction) = {1}; index_end(direction) = {dim(direction)};
- 将
index_first
在指定direction
方向上的索引设置为1
,表示**只选取该维度上的第一个元素。**例如,如果direction = 2
且X
是一个 3D 张量,index_first
将变为{'', '1', ''}
,表示选取第二维的第一个元素。 - 将
index_end
在指定direction
方向上的索引设置为该维度的最后一个元素,dim(direction)
返回张量X
在direction
****方向上的尺寸。例如,如果direction = 2
且dim = [4, 5, 6]
,index_end
将变为{'', '5', ''}
,表示选取第二维的最后一个元素。
- 将
-
计算起始和结束切片的差
slice = X(index_first{:}) - X(index_end{:});
- 计算
X
在direction
方向上第一个元素与最后一个元素的差,形成一个切片slice
。
- 计算
-
计算沿指定方向的差分
DX = diff(X, 1, direction);
- 使用 MATLAB 的
diff
函数计算沿direction
方向的一阶差分,结果存储在DX
中。
- 使用 MATLAB 的
-
将差分结果和切片差拼接起来
DX = cat(direction, DX, slice);
- 将差分结果
DX
与切片差slice
沿direction
方向拼接,得到最终的差分张量DX
。
- 将差分结果
通过这些步骤,
porder_diff
函数实现了沿指定方向计算张量的差分,并将起始和结束元素的差作为补充,确保差分张量的维度与原始张量一致。 -
porder_diff_T
- 计算张量
Y
在指定方向index
****上的转置差分。 - 应返回一个与
Y
尺寸相同的张量,表示在该方向上的转置差分结果。
diff_element
- 计算在指定方向
direction
上的差分元素矩阵。 - 应返回一个与
dim
兼容的张量,表示在该方向上的差分元素。
function Eny = diff_element(dim, direction)
d = length(dim); % 获取维度向量的长度
e = ones(1, d); % 创建一个长度为 d 的全为 1 的向量
element1 = ones(e); % 创建一个形状为 e 的全 1 数组
element2 = -1 * ones(e); % 创建一个形状为 e 的全 -1 数组
element = cat(direction, element1, element2); % 在指定方向上拼接 element1 和 element2
Eny = (abs(psf2otf(element, dim))).^2; % 将点扩散函数转换为光学传递函数并取绝对值平方
end
prox_htnn_F
- 计算张量
X
在傅里叶变换下的近端算子,参数tau
用于调整。 - 应返回一个与
X
尺寸相同的张量,表示近端算子的结果和张量核范数。
function [X,htnn,tsvd_rank] = prox_htnn_F(Y,rho)
%The proximal operator for the order-D tensor nuclear norm under Discrete Fourier Transform (DFT)
p = length(size(Y));
n = zeros(1,p);
for i = 1:p
n(i) = size(Y,i);
end
X = zeros(n);
L = ones(1,p);
for i = 3:p
Y = fft(Y,[],i);
L(i) = L(i-1) * n(i);
end
htnn = 0;
tsvd_rank = 0;
[U,S,V] = svd(Y(:,:,1),'econ');
S = diag(S);
r = length(find(S>rho));
if r>=1
S = max(S(1:r)-rho,0);
X(:,:,1) = U(:,1:r)*diag(S)*V(:,1:r)';
htnn = htnn+sum(S);
tsvd_rank = max(tsvd_rank,r);
end
for j = 3 : p
for i = L(j-1)+1 : L(j)
%
I = unfoldi(i,j,L);
halfnj = floor(n(j)/2)+1;
%
if I(j) <= halfnj && I(j) >= 2
[U,S,V] = svd(Y(:,:,i),'econ');
S = diag(S);
r = length(find(S>rho));
if r>=1
S = max(S(1:r)-rho,0);
X(:,:,i) = U(:,1:r)*diag(S)*V(:,1:r)';
htnn = htnn+sum(S)*2;
tsvd_rank = max(tsvd_rank,r);
end
%Conjugation property
elseif I(j) > halfnj
%
n_ = nc(I,j,n);
%
i_ = foldi(n_,j,L);
X(:,:,i) = conj( X(:,:,i_));
end
end
end
htnn = htnn/prod(n(3:end));
for i = p:-1:3
X = (ifft(X,[],i));
end
X = real(X);
prox_htnn_C
- 计算张量
X
在余弦变换下的近端算子,参数tau
用于调整。 - 应返回一个与
X
尺寸相同的张量,表示近端算子的结果和张量核范数。
⭐TCTV_TC
默认参数设置
dim = size(M); % 获取张量M的尺寸,存储在dim中。
d = ndims(M); % 获取张量M的阶数(维度数),存储在d中。
transform = 'DFT'; % 设置默认的变换类型为DFT(离散傅里叶变换)。
% 对于从第三维度到第d维度的每一维度,生成对应尺寸的DFT矩阵,并存储在transform_matrices中。
transform_matrices = arrayfun(@(x) dftmtx(dim(x)), 3:d, 'UniformOutput', false);
directions = 1:d; % 设置默认的方向为从1到d的所有维度。
tol = 1e-8; % 设置默认的终止容差为1e-8。
max_iter = 500; % 设置默认的最大迭代次数为500次。
rho = 1.1; % 设置默认的步长增加比例rho为1.1。
mu = 1e-4; % 设置默认的初始步长mu为1e-4。
max_mu = 1e10; % 设置默认的最大步长max_mu为1e10。
detail = 1; % 设置默认的详细信息显示标志为1(显示详细信息)。
- 详细解释:
dim = size(M);
:这个语句获取张量 M 的尺寸。例如,如果M是一个3维张量,其尺寸为[3, 4, 5],那么 dim 将是 [3, 4, 5]。d = ndims(M);
:这个语句获取张量M的维数(阶数)。例如,如果M是一个3维张量,d将是3。transform = 'DFT';
:设定默认的变换类型为DFT,表示在高阶TSVD中使用离散傅里叶变换。transform_matrices = arrayfun(@(x) dftmtx(dim(x)), 3:d, 'UniformOutput', false);
:- 这里使用了
arrayfun
函数对从第三维度到第 d 维度的每一个维度生成对应的DFT矩阵。 @(x) dftmtx(dim(x))
:这是一个匿名函数,生成一个尺寸为dim(x)
的DFT矩阵。'UniformOutput', false
:这是一个参数,指示arrayfun
返回的结果应该是一个单元数组(cell array),因为每个函数的输出可能是不同的大小或类型。
- 这里使用了
使用用户选项覆盖默认设置
if exist('opts', 'var')
if isfield(opts, 'transform'), transform = opts.transform; end
if isfield(opts, 'transform_matrices'), transform_matrices = opts.transform_matrices; end
if isfield(opts, 'directions'), directions = opts.directions; end
if isfield(opts, 'tol'), tol = opts.tol; end
if isfield(opts, 'max_iter'), max_iter = opts.max_iter; end
if isfield(opts, 'rho'), rho = opts.rho; end
if isfield(opts, 'mu'), mu = opts.mu; end
if isfield(opts, 'max_mu'), max_mu = opts.max_mu; end
if isfield(opts, 'detail'), detail = opts.detail; end
end
对变量初始化
% 变量初始化
n = length(directions); % 获取**方向的数量**,即考虑的局部平滑方向的数量。
X = randn(dim); % 初始化 X 为与 M 相同尺寸的随机张量。
X(Omega) = M(Omega); % 将观测值的索引处的 X 的值设置为 M 的观测值。
E = zeros(dim); % 初始化E为与M相同尺寸的零张量。
Lambda = zeros(dim); % 初始化Lambda为与M相同尺寸的零张量。
G = cell(1, n); % 初始化G为包含n个元素的单元数组,用于**存储各方向上的差分结果。**
Gamma = cell(1, n); % 初始化Gamma为包含n个元素的单元数组,用于**存储各方向上的对偶变量。**
% 逐方向初始化G和Gamma
for i = 1:n
index = directions(i); % 获取当前方向的索引。
G{index} = porder_diff(X, index); % 计算X在当前方向上的差分,存储在G中。
Gamma{index} = zeros(dim); % 初始化当前方向上的 Gamma 为与M相同尺寸的零张量。
end
- 详细解释:
n = length(directions);
- 获取方向的数量,即考虑的局部平滑方向的数量。
directions
包含了所有要考虑的方向的索引[1:d]
,这里n
就是这些索引的数量。
- 获取方向的数量,即考虑的局部平滑方向的数量。
X = randn(dim);
- 初始化一个与
M
尺寸相同的随机张量X
,元素服从标准正态分布。
- 初始化一个与
E = zeros(dim);
- 初始化与
M
相同尺寸的零张量E
。E
用来存储误差。
- 初始化与
Lambda = zeros(dim);
- 初始化与
M
相同尺寸的零张量Lambda
。Lambda
是 ADMM 算法中的对偶变量。
- 初始化与
G = cell(1, n);
- 初始化包含
n
个元素的单元数组G
,用来存储各方向上的差分结果。
- 初始化包含
Gamma = cell(1, n);
- 初始化包含
n
个元素的单元数组Gamma
,用来存储各方向上的对偶变量。
- 初始化包含
G{index} = porder_diff(X, index);
- 计算张量
X
在当前方向上的差分,并存储在G
中。porder_diff
是一个函数,用于计算指定方向上的差分。
- 计算张量
Gamma{index} = zeros(dim);
- 初始化当前方向上的对偶变量
Gamma
为与M
相同尺寸的零张量。
- 初始化当前方向上的对偶变量
FFT 设置
% FFT 设置
T = zeros(dim); % 初始化与 M 尺寸相同的零张量 T。
for i = 1:n
Eny = diff_element(dim, directions(i)); % 计算**在当前方向上的差分元素矩阵 Eny。**
T = T + Eny; % 将 Eny 累加到 T 中。
end
-
详细解释:
T = zeros(dim);
- 初始化与
M
尺寸相同的零张量T
。这个张量将用于存储累加的差分元素矩阵。
- 初始化与
for i = 1:n
- 循环遍历每个方向的索引。
Eny = diff_element(dim, directions(i));
- 计算在当前方向上的差分元素矩阵
Eny
。diff_element
是一个函数,用于计算指定方向上的差分元素矩阵。 dim
是张量M
的尺寸。directions(i)
是当前方向的索引。
- 计算在当前方向上的差分元素矩阵
T = T + Eny;
- 将当前方向的差分元素矩阵
Eny
累加到T
中。
- 将当前方向的差分元素矩阵
通过这些步骤,
T
将包含所有指定方向上的差分元素矩阵的累加结果。这个累加的差分元素矩阵在后续的 FFT 运算中将会被用到,以调整X
的更新过程。
主循环
iter = 0;
while iter < max_iter
iter = iter + 1;
Xk = X;
Ek = E;
% 更新 X
H = zeros(dim);
for i = 1:n
index = directions(i);
H = H + porder_diff_T(mu * G{index} - Gamma{index}, index);
end
X = real(ifftn(fftn(mu * (M - E) + Lambda + H) ./ (mu * (1 + T))));
% 更新 G
for i = 1:n
index = directions(i);
switch transform
case 'DFT'
[G{index}, tnn_G{index}] = prox_htnn_F(porder_diff(X, index) + Gamma{index} / mu, 1 / (n * mu));
case 'DCT'
[G{index}, tnn_G{index}] = prox_htnn_C(porder_diff(X, index) + Gamma{index} / mu, 1 / (n * mu));
otherwise % 自定义变换
[G{index}, tnn_G{index}] = prox_htnn_C(transform_matrices, porder_diff(X, index) + Gamma{index} / mu, 1 / (n * mu));
end
end
% 更新 E
E = M - X + Lambda / mu;
E(Omega) = 0;
% 停止准则
dY = M - X - E;
chgX = max(abs(Xk(:) - X(:)));
chgE = max(abs(Ek(:) - E(:)));
chg = max([chgX, chgE, max(abs(dY(:)))]);
if chg < tol, break; end
% 显示更新细节
if detail && (iter == 1 || mod(iter, 10) == 0)
obj = sum(cell2mat(tnn_G)) / n;
err = norm(dY(:), 'fro');
disp(['iter ', num2str(iter), ', mu=', num2str(mu), ', obj=', num2str(obj), ', err=', num2str(err)]);
end
% 更新乘子和 mu
Lambda = Lambda + mu * dY;
for i = 1:n
index = directions(i);
Gamma{index} = Gamma{index} + mu * (porder_diff(X, index) - G{index});
end
mu = min(rho * mu, max_mu);
end
end