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
问题描述
我们希望通过最小化张量的相关全变差（Tensor Correlated Total Variation, TCTV）范数来完成张量的恢复。这个问题可以用如下优化模型表示：
min X ∥ X ∥ TCTV subject to P Ω ( X ) = P Ω ( M ) \min_{X} \|X\|_{\text{TCTV}} \quad \text{subject to} \quad P_\Omega(X) = P_\Omega(M) Xmin∥X∥TCTVsubject toPΩ(X)=PΩ(M)
其中：
- M M M 是观测到的 p p p 阶张量。
- X X X 是恢复的 p p p 阶张量。
- ∥ X ∥ TCTV \|X\|_{\text{TCTV}} ∥X∥TCTV 是张量 X X X 的 TCTV 范数。
- P Ω ( ⋅ ) P_\Omega(\cdot) PΩ(⋅) 是在已知观测位置集 Ω \Omega Ω 上的投影操作。
更新公式的推导
1. 优化 X X X
在更新 X X X 时，我们需要最小化如下目标函数：
∥ X ∥ TCTV + μ 2 ∥ P Ω ( X ) − P Ω ( M ) ∥ F 2 + μ 2 ∥ M − X − E + Λ μ ∥ F 2 \|X\|_{\text{TCTV}} + \frac{\mu}{2} \|P_\Omega(X) - P_\Omega(M)\|_F^2 + \frac{\mu}{2} \|M - X - E + \frac{\Lambda}{\mu}\|_F^2 ∥X∥TCTV+2μ∥PΩ(X)−PΩ(M)∥F2+2μ∥M−X−E+μΛ∥F2
为了使这部分容易处理，使用频域变换（例如傅里叶变换）将问题转换到频域。通过这种方式，差分操作可以转化为简单的乘法，从而简化了 TCTV范数的计算。
X = F − 1 ( F ( μ ( M − E + Λ / μ ) + H ) μ ( 1 + T ) ) X = \mathcal{F}^{-1} \left( \frac{\mathcal{F}(\mu (M - E + \Lambda / \mu) + H)}{\mu(1 + T)} \right) X=F−1(μ(1+T)F(μ(M−E+Λ/μ)+H))
其中 F \mathcal{F} F 表示傅里叶变换， F − 1 \mathcal{F}^{-1} F−1 表示逆傅里叶变换， T T T 是差分操作频谱。
2. 优化 E
在更新 E E E 时，我们需要最小化如下目标函数：
E = arg min E ( μ 2 ∥ E − ( M − X + Λ μ ) ∥ F 2 ) E = \arg\min_E \left( \frac{\mu}{2} \|E - (M - X + \frac{\Lambda}{\mu})\|_F^2 \right) E=argEmin(2μ∥E−(M−X+μΛ)∥F2)
这是一个简单的二次优化问题，其解析解为：
E = M − X + Λ / μ E = M - X + \Lambda / \mu E=M−X+Λ/μ
为了确保
E
E
E 在
Ω
\Omega
Ω 上等于零：
E
(
Ω
)
=
0
E(\Omega) = 0
E(Ω)=0
3. 更新拉格朗日乘子 Λ \Lambda Λ
更新拉格朗日乘子的公式源自标准的 ADMM 更新规则：
Λ
k
+
1
=
Λ
k
+
μ
(
M
−
X
k
+
1
−
E
k
+
1
)
\Lambda^{k+1} = \Lambda^k + \mu (M - X^{k+1} - E^{k+1})
Λk+1=Λk+μ(M−Xk+1−Ek+1)
具体实现
这些公式在代码中通过以下步骤实现：
-
变量初始化:
n = length(directions); X = randn(dim); X(Omega) = M(Omega); E = zeros(dim); Lambda = zeros(dim);- X X X : 初始恢复张量，通常初始化为在 Ω \Omega Ω 上等于 M M M 的随机张量。
- E E E : 初始误差张量，初始化为零张量。
- Λ \Lambda Λ : 拉格朗日乘子，初始化为零张量。
- G G G : 每个方向的差分张量，初始化为零张量。
- Γ \Gamma Γ : 每个方向的拉格朗日乘子，初始化为零张量。s
- 其他参数：根据输入
opts初始化。
-
更新 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))));X k + 1 = arg min X ( ∥ X ∥ TCTV + μ 2 ∥ P Ω ( X ) − P Ω ( M ) ∥ F 2 + μ 2 ∥ M − X − E + Λ μ ∥ F 2 ) X^{k+1} = \arg\min_X \left( \|X\|_{\text{TCTV}} + \frac{\mu}{2} \|P_\Omega(X) - P_\Omega(M)\|_F^2 + \frac{\mu}{2} \|M - X - E + \frac{\Lambda}{\mu}\|_F^2 \right) Xk+1=argXmin(∥X∥TCTV+2μ∥PΩ(X)−PΩ(M)∥F2+2μ∥M−X−E+μΛ∥F2)
-
更新 E:
E = M - X + Lambda / mu; E(Omega) = 0;E k + 1 = arg min E ( μ 2 ∥ E − ( M − X + Λ μ ) ∥ F 2 ) E^{k+1} = \arg\min_E \left( \frac{\mu}{2} \|E - (M - X + \frac{\Lambda}{\mu})\|_F^2 \right) Ek+1=argEmin(2μ∥E−(M−X+μΛ)∥F2)
-
更新拉格朗日乘子 Λ \Lambda Λ:
Lambda = Lambda + mu * (M - X - E); 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);Λ k + 1 = Λ k + μ ( M − X k + 1 − E k + 1 ) \Lambda^{k+1} = \Lambda^k + \mu (M - X^{k+1} - E^{k+1}) Λk+1=Λk+μ(M−Xk+1−Ek+1)
⭐TCTV_TRPCA
问题描述
我们希望通过最小化张量的相关全变差（Tensor Correlated Total Variation, TCTV）范数和稀疏误差项来从观察到的张量 M M M 中分离出低秩张量 X X X 和稀疏张量 E E E 。这个问题可以用如下优化模型表示：
min X , E ∥ X ∥ TCTV + λ ∥ E ∥ 1 subject to M = X + E \min_{X, E} \|X\|_{\text{TCTV}} + \lambda \|E\|_1 \quad \text{subject to} \quad M = X + E X,Emin∥X∥TCTV+λ∥E∥1subject toM=X+E
其中：
- M M M 是观测到的 p p p 阶张量。
- X X X 是恢复的低秩 p p p 阶张量。
- E E E 是稀疏 p p p 阶张量。
- ∥ X ∥ TCTV \|X\|_{\text{TCTV}} ∥X∥TCTV 是张量 X X X 的 TCTV 范数。
- ∥ E ∥ 1 \|E\|_1 ∥E∥1 是张量 E E E 的 ℓ 1 \ell_1 ℓ1 范数。
- λ \lambda λ 是控制低秩项和稀疏项之间权衡的正则化参数。
相关公式
-
张量的相关全变差 (TCTV) 范数:
∥ X ∥ TCTV = ∑ i = 1 n ∥ ∇ i X ∥ ∗ \|X\|_{\text{TCTV}} = \sum_{i=1}^n \| \nabla_i X \|_* ∥X∥TCTV=i=1∑n∥∇iX∥∗
其中 ∇ i \nabla_i ∇i 是在第 i i i 个方向上的差分操作， ∥ ⋅ ∥ ∗ \| \cdot \|_* ∥⋅∥∗ 表示张量的核范数。
-
稀疏张量的 ℓ 1 \ell_1 ℓ1 范数:
∥ E ∥ 1 = ∑ i = 1 p ∑ j = 1 n i ∣ E i j ∣ \|E\|_1 = \sum_{i=1}^p \sum_{j=1}^{n_i} |E_{ij}| ∥E∥1=i=1∑pj=1∑ni∣Eij∣
各部分代码的数学模型对应关系
-
变量初始化:
X = randn(dim); E = zeros(dim); Lambda = zeros(dim);- X X X : 初始低秩张量，通常初始化为随机张量。
- E E E : 初始稀疏张量，初始化为零张量。
- Λ \Lambda Λ : 拉格朗日乘子，初始化为零张量。
- G G G : 每个方向的差分张量，初始化为零张量。
- Γ \Gamma Γ : 每个方向的拉格朗日乘子，初始化为零张量。
- 其他参数：根据输入
opts初始化。
-
更新 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))));X = F − 1 ( F ( μ ( M − E + Λ / μ ) + H ) μ ( 1 + T ) ) X = \mathcal{F}^{-1} \left( \frac{\mathcal{F}(\mu(M - E + \Lambda/\mu) + H)}{\mu(1 + T)} \right) X=F−1(μ(1+T)F(μ(M−E+Λ/μ)+H))
其中 F \mathcal{F} F 表示快速傅里叶变换， H H H 是由 G G G 和 Γ \Gamma Γ 计算的辅助张量， T T 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)); case 'other' [G{index}, tnn_G{index}] = prox_htnn_C(transform_matrices, porder_diff(X, index) + Gamma{index} / mu, 1 / (n * mu)); end endG i = prox htnn ( ∇ i X + Γ i / μ , 1 n μ ) G_i = \text{prox}_{\text{htnn}}\left( \nabla_i X + \Gamma_i / \mu, \frac{1}{n\mu} \right) Gi=proxhtnn(∇iX+Γi/μ,nμ1)
其中 ∇ i \nabla_i ∇i 表示在第 i i i 个方向上的差分操作。
-
更新 E:
E = prox_l1(M - X + Lambda / mu, lambda / mu);E = prox ℓ 1 ( M − X + Λ / μ , λ / μ ) E = \text{prox}_{\ell_1}(M - X + \Lambda / \mu, \lambda / \mu) E=proxℓ1(M−X+Λ/μ,λ/μ)
-
**停止条件：**当迭代步数达到最大迭代次数或变化量小于预设容差时停止。
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 -
更新拉格朗日乘子和步长:
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);Λ = Λ + μ ( M − X − E ) \Lambda = \Lambda + \mu (M - X - E) Λ=Λ+μ(M−X−E)
Γ i = Γ i + μ ( ∇ i X − G i ) \Gamma_i = \Gamma_i + \mu (\nabla_i X - G_i) Γi=Γi+μ(∇iX−Gi)
μ = min ( ρ μ , max_mu ) \mu = \min(\rho \mu, \text{max\_mu}) μ=min(ρμ,max_mu)
4876
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
