【张量｜TCTV】t-CTV 代码分析与实验结果

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);

详细解释：
1. 获取张量 X 的尺寸
```
dim = size(X);
```
  - dim 是一个包含 X 各个维度尺寸的向量。例如，如果 X 是一个 3D 张量，其尺寸为 [4, 5, 6]，那么 dim 将是 [4, 5, 6]。
2. 初始化索引
```
index_first = repmat({':'}, 1, ndims(X));
index_end = repmat({':'}, 1, ndims(X));
```
  - index_first 和 index_end 初始化为适用于多维数组的全选索引。例如，对于一个 3D 张量，这些索引初始为 {':', ':', ':'}。表示全选该维度上的所有元素。
3. 设置起始和结束索引
```
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', ''}，表示选取第二维的最后一个元素。
4. 计算起始和结束切片的差
```
slice = X(index_first{:}) - X(index_end{:});
```
  - 计算 X 在 direction 方向上第一个元素与最后一个元素的差，形成一个切片 slice。
5. 计算沿指定方向的差分
```
DX = diff(X, 1, direction);
```
  - 使用 MATLAB 的 diff 函数计算沿 direction 方向的一阶差分，结果存储在 DX 中。
6. 将差分结果和切片差拼接起来
```
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\|_{\text{TCTV}} \quad \text{subject to} \quad P_\Omega(X) = P_\Omega(M)$

其中：

$M$ 是观测到的 $p$ 阶张量。
$X$ 是恢复的 $p$ 阶张量。
$\|X\|_{\text{TCTV}}$ 是张量 $X$ 的 TCTV 范数。
$P_\Omega(\cdot)$ 是在已知观测位置集 $\Omega$ 上的投影操作。

更新公式的推导

1. 优化 $X$

在更新 $X$ 时，我们需要最小化如下目标函数：

$\|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$

为了使这部分容易处理，使用频域变换（例如傅里叶变换）将问题转换到频域。通过这种方式，差分操作可以转化为简单的乘法，从而简化了 TCTV范数的计算。

$\mathcal{F}^{-1} \left( \frac{\mathcal{F}(\mu (M - E + \Lambda / \mu) + H)}{\mu(1 + T)} \right)$

其中 $\mathcal{F}$ 表示傅里叶变换， $\mathcal{F}^{-1}$ 表示逆傅里叶变换， $T$ 是差分操作频谱。

2. 优化 E

在更新 $E$ 时，我们需要最小化如下目标函数：

$\arg\min_E \left( \frac{\mu}{2} \|E - (M - X + \frac{\Lambda}{\mu})\|_F^2 \right)$

这是一个简单的二次优化问题，其解析解为：

$\Lambda / \mu$

为了确保 $E$ 在 $\Omega$ 上等于零：
$E(\Omega) = 0$

3. 更新拉格朗日乘子 $\Lambda$

更新拉格朗日乘子的公式源自标准的 ADMM 更新规则：
$\Lambda^{k+1} = \Lambda^k + \mu (M - X^{k+1} - E^{k+1})$

具体实现

这些公式在代码中通过以下步骤实现：

变量初始化:
```
n = length(directions);
X = randn(dim);
X(Omega) = M(Omega);
E = zeros(dim);
Lambda = zeros(dim);
```
- $X$ : 初始恢复张量，通常初始化为在 $\Omega$ 上等于 $M$ 的随机张量。
- $E$ : 初始误差张量，初始化为零张量。
- $\Lambda$ : 拉格朗日乘子，初始化为零张量。
- $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 \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)$
更新 E:
```
E = M - X + Lambda / mu;
E(Omega) = 0;
```
$E^{k+1} = \arg\min_E \left( \frac{\mu}{2} \|E - (M - X + \frac{\Lambda}{\mu})\|_F^2 \right)$
更新拉格朗日乘子 $\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);
```
$\Lambda^{k+1} = \Lambda^k + \mu (M - X^{k+1} - E^{k+1})$

⭐TCTV_TRPCA

问题描述

我们希望通过最小化张量的相关全变差（Tensor Correlated Total Variation, TCTV）范数和稀疏误差项来从观察到的张量 $M$ 中分离出低秩张量 $X$ 和稀疏张量 $E$ 。这个问题可以用如下优化模型表示：

$\min_{X, E} \|X\|_{\text{TCTV}} + \lambda \|E\|_1 \quad \text{subject to} \quad M = X + E$

其中：

$M$ 是观测到的 $p$ 阶张量。
$X$ 是恢复的低秩 $p$ 阶张量。
$E$ 是稀疏 $p$ 阶张量。
$\|X\|_{\text{TCTV}}$ 是张量 $X$ 的 TCTV 范数。
$E\|_1$ 是张量 $E$ 的 $\ell_1$ 范数。
$\lambda$ 是控制低秩项和稀疏项之间权衡的正则化参数。

各部分代码的数学模型对应关系

变量初始化:
```
X = randn(dim);
E = zeros(dim);
Lambda = zeros(dim);
```
- $X$ : 初始低秩张量，通常初始化为随机张量。
- $E$ : 初始稀疏张量，初始化为零张量。
- $\Lambda$ : 拉格朗日乘子，初始化为零张量。
- $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))));
```
$\mathcal{F}^{-1} \left( \frac{\mathcal{F}(\mu(M - E + \Lambda/\mu) + H)}{\mu(1 + T)} \right)$

其中 $\mathcal{F}$ 表示快速傅里叶变换， $H$ 是由 $G$ 和 $\Gamma$ 计算的辅助张量， $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
end
```
$G_i = \text{prox}_{\text{htnn}}\left( \nabla_i X + \Gamma_i / \mu, \frac{1}{n\mu} \right)$

其中 $\nabla_i$ 表示在第 $i$ 个方向上的差分操作。
更新 E:
```
E = prox_l1(M - X + Lambda / mu, lambda / mu);
```
$\text{prox}_{\ell_1}(M - X + \Lambda / \mu, \lambda / \mu)$

**停止条件：**当迭代步数达到最大迭代次数或变化量小于预设容差时停止。

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);
```
$\Lambda = \Lambda + \mu (M - X - E)$

$\Gamma_i = \Gamma_i + \mu (\nabla_i X - G_i)$

$max_mu ) \mu = \min(\rho \mu, \text{max\_mu})$