Complexity of solving an image differential linear system

Define an "image differential linear system" as a linear system $$A\mathbf{x}=\mathbf{b}$$ wherein $$\mathbf{x}$$ contains the ($$\mathbb{R}$$) pixels of an image and each row of $$A$$ constrains the difference between pixels that neighbor each other (vertically or horizontally) in the image.

Such a linear system:

• has entries $$A_{ij}\in\{-1,0,1\} \forall i,j$$
• is very sparse (having exactly one $$-1$$ and one $$+1$$ entry in each row
• is overdetermined (having almost twice as many rows as columns)
• is "planar" in the sense that its rows can be seen as edges in a planar graph (grid in fact) of pixels

My question is: Can we say anything about the complexity of solving such a linear system directly?

• By "solving" I mean finding the least-squares solution (as the system is over-determined).
• By "directly" I mean not using iterated approximations/relaxations.

As practical images normally have millions of pixels I am looking to gain whatever intuitions/insights I can on the asymptotic complexity of solving this commonly-occurring linear system.

• Solving it with which algorithm? Dec 1, 2020 at 18:26
• @FedericoPoloni Whichever algorithm would be the fastest for this type of sparse system. (Multigrid relaxation, with supposedly linear complexity, is often used but that is iterative method, not direct.) Dec 1, 2020 at 19:30
• Is your question "which is the best algorithm to solve these systems" then, essentially? Asking for the complexity seems a red herring. Dec 1, 2020 at 20:06
• @FedericoPoloni Yes nearly. I'm asking what is a achievable with direct methods, with a view to revisiting a decision I had made some time ago to use (supposedly O(n)) iterative approximation/relaxation. Dec 1, 2020 at 20:14
• This is about 2D graphs, not images only. I would rewrite in terms of graph theory. Note that $A$ is an incidence matrix. Take a look at this. Dec 2, 2020 at 19:14

TL;DR

For an image of size $$m\times n$$ you can solve this problem in $$O(nm(\log(n) + \log(m)))$$. In fact, there is nothing to "solve", the solution can be written down analytically. The differencing operator is linear space invariant, i.e. a convolution, and you are asking how to do a "de-convolution". This can be done efficiently in the frequency domain: transform the image with a 2D FFT and divide by the transfer function of the differencing filter, then transform back. There are some details that need to be addressed: the convolution is not cyclic, and the solution you ask for is a least-squares one. This is explained below, where I show the relation to finite-difference approximations of the Poisson equation and Sylvester equations.

Details

To make the notations simpler, let's assume $$m = n$$. Using the finite-difference matrix $$D = \begin{bmatrix}1 &-1& & & \\ & 1 & -1 & & \\ & & \ddots & \ddots & \\ & & & 1 & -1\end{bmatrix},$$ we can write your equation as $$\begin{bmatrix}I \otimes D\\ D^T \otimes I\end{bmatrix} x\approx \begin{bmatrix}b_r\\b_c\end{bmatrix},$$ which is to be solved in the least-squares sense, with $$b_r, b_c$$ the differences in the row and column directions. The normal equations for the least-squares problem are obtained by multiplying both sides by $$\left[I\otimes D^T, \ D\otimes I\right]$$, which gives $$(I\otimes D^TD + DD^T\otimes I)x=(I\otimes D^T) b_r + (D\otimes I)b_c.$$ This can be written as a Sylvester equation (see https://en.wikipedia.org/wiki/Kronecker_product), $$A X + XB = C,$$ where $$A = D^TD$$, $$B = DD^T$$ and $$C = DB_r + B_cD$$, and the upper-case letters are matrix versions of the lower-case letters representing column vectors. In particular, $$X$$ is just the original $$n\times n$$ image. We note the $$A$$ and $$B$$ are finite-difference 1D Laplacian matrices, so your problem is reduced to a discrete version of the Poisson equation.

In your question you did not mention any boundary conditions, without which the problem is not well defined. If you only specify the differences between neighboring pixels, you can add an arbitrary constant to your solution without affecting the least-squares error. There are several possible boundary-conditions to choose from, for example, you could use Dirichlet b.c. for all edges. Choosing Neumann b.c. for all edges will not yield a unique solution, so that should be avoided. Whichever b.c. you choose, it will determine the exact form of the $$A, \ B$$ (and $$C$$) matrices, which will be slightly different from $$DD^T$$ and $$D^TD$$. For the solution to be unique the spectrum of $$A$$ and $$-B$$ must be disjoint (https://en.wikipedia.org/wiki/Sylvester_equation).

For all common choices of b.c., the resulting $$A$$ and $$B$$ have known eigenvalues and eigenvectors (https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors_of_the_second_derivative), and we may write, $$A = U_A\Lambda_A U_A^T, \ B = U_B\Lambda_B U_B^T,$$ where $$\Lambda_A = \mathrm{diag}(\lambda_A), \Lambda_B = \mathrm{diag}(\lambda_B)$$, and $$U_A^T U_A = I, U_B^T U_B = I$$. Moreover, every choice of b. c. is associated with some version of the discrete cosine or sine transforms for which fast algorithms exist, so that computing the product of $$U_A, \ U_B$$ times a vector can be done in $$O(n \log n)$$ (https://en.wikipedia.org/wiki/Discrete_sine_transform).

Inserting the diagonalizations into the Sylvester equation we have, $$U_A\Lambda_A U_A^T X + X U_B\Lambda_B U_B^T = C \implies \Lambda_A (U_A^T X U_B) + (U_A^TX U_B)\Lambda_B = U_A^T C U_B,$$ which is a Sylvester equation in the transformed image $$\hat{X} = U_A^T X U_B$$ with diagonal coefficient matrices. Such equations have a known solution, $$\hat{X}_{ij} = \frac{(U_A^T C U_B)_{ij}}{{(\lambda_A)}_i + {(\lambda_B)}_j}$$ (note the relation to the condition for the existence of the solution to the Sylvester equation). Lastly, transform back to the image to obtain, $$X = U_A \hat{X} U_B^T.$$