1
$\begingroup$

I'm reading a paper - Total Variation Superiorized Conjugate Gradient Method for Image Reconstruction on total variation regularization and conjugate gradients. In page $3$, the authors define the total variation norm of a (vectorized) image as the $l_1$ norm applied to a matrix-vector product. More precisely, they define the function \begin{align} TV(x)=\|Rx\|_1, \end{align} where the input image $x \in \mathbb{R}^n$ is vectorized. In practice, I know how to calculate $TV(x)$ procedurally using \begin{align} \sum_{i,j} \left| X_{i,j} - X_{i-1,j} \right| + \left| X_{i,j} - X_{i,j-1} \right| \end{align} where $X=\text{devec}(x)$, and which I believe is the same thing as the authors' $TV(x)$ (right?). Later in the paper, the authors approximate $\nabla TV (x)$ using an expression that involves the matrix $R$ and its transpose $R^T$.

So my question is: How can I construct such a matrix $R$ so that I can transpose it? or even better: can I calculate the matrix-vector products $Rx$ and $R^Tx$ without having to construct the matrix itself?

$\endgroup$

1 Answer 1

1
$\begingroup$

I just found out on my own. Assuming that $X \in \mathbb{R}^{m \times n}$, we define the difference matrices $\Delta_I X = [ X_{i,j} - X_{i-1,j} ]$ and $\Delta_J X = [X_{i,j} - X_{i,j-1}]$, where negative indices are defined $\text{mod}\ m$ or $\text{mod}\ n$ for $i$ or $j$, respectivelly. I first noticed by inspection that both $\Delta_I X$ and $\Delta_J X$ can be obtained via the multiplication of $X$ with circulant matrices $R_I$ and $R_J$ given by \begin{align} R_I &= \begin{pmatrix} 1 & 0 & 0 &\dots & 0 & -1 \\ -1 & 1 & 0 & \dots & 0 & 0 \\ 0 & -1 & 1 & \dots & 0 & 0 \\ & & & \vdots & & \\ 0 & 0 & 0 & \dots & -1 & 1 \end{pmatrix} \in \mathbb{R}^{m \times m} \\ \end{align} and \begin{align} R_J &= \begin{pmatrix} 1 & -1 & 0 &\dots & 0 \\ 0 & 1 & -1 & \dots & 0 \\ 0 & 0 & 1 & \dots & 0 \\ & & & \vdots & \\ 0 & 0 & 0 & \dots & -1 \\ -1 & 0 & 0 & \dots & 1 \end{pmatrix} \in \mathbb{R}^{n \times n}. \end{align} Precisely, we have that $\Delta_I X = R_I X$ and $\Delta_J X = X R_J$. Using the Kronecker Product $\otimes$ we can write these expressions in vector form as \begin{align} \text{vec}(\Delta_I X) &= (R_I \otimes I_n) \text{vec}(X) = (R_I \otimes I_n) x\\ \text{vec}(\Delta_J X) &= (I_m \otimes R_J) \text{vec}(X) = (I_m \otimes R_J) x. \end{align} Concatenating $R_I \otimes I_n$ and $I_m \otimes R_J$ we have the desired matrix \begin{align} R &= \begin{pmatrix} R_I \otimes I_n \\ I_m \otimes R_J \end{pmatrix} \in \mathbb{R}^{2mn \times mn}. \end{align} Obviously, the matrix-vector product $Rx$ can be implemented indirectly, that is, without having to store $R$ in memory. This can be done by shifting and subtracting. In any case, this matrix representation of the "pre TV function" is very interesting from an analytical standpoint. Since transposition distributes over the Kronecker product, we have that \begin{align} R^T &= \begin{pmatrix} R_I^T \otimes I_n & I_m \otimes R_J^T \end{pmatrix} \in \mathbb{R}^{mn \times 2mn}. \end{align} I also noticed by inspection that \begin{align} \text{devec}[(R_I^T \otimes I_n)x] &= R_I^T X = \text{flip}_J \left( R_I \text{flip}_J (X) \right), \tag{$1$} \end{align} where $\text{flip}_J(M)$ is the matrix $M$ with all columns written in reversed order. Similarly, we have that \begin{align} \text{devec}[(I_m \otimes R_J^T)x] &= X R_J^T = \text{flip}_I \left(\text{flip}_I(X) R_J \right), \tag{$2$} \end{align} where $\text{flip}_J(M)$ is the matrix $M$ with all rows written in reversed order. Therefore, using $(1)$ and $(2)$, the matrix-vector product $R^Tx$ can be easily implemented with flips, shifts and subtractions similar to $Rx$, and without the need of storing $R$ nor $R^T$ in memory.

$\endgroup$
2
  • 2
    $\begingroup$ Pay attention to different variants of the TV Norm for 2D signals: dsp.stackexchange.com/questions/57977. $\endgroup$
    – Royi
    Apr 12, 2023 at 7:59
  • $\begingroup$ Thank you very much! $\endgroup$
    – mlbj
    Apr 12, 2023 at 13:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.