Computing adjoint mismatch of linear maps

This paper proposes and analyzes a stochastic algorithm that provably converges almost surely to the operator norm of the difference between two linear maps—one accessible only via forward evaluation and the other only via adjoint evaluation—by performing a random search on a generalized Rayleigh quotient with optimal step sizes.

The Gist

Here is an explanation of the paper "Computing adjoint mismatch of linear maps," translated into simple language with creative analogies.

The Big Picture: The "Black Box" Problem

Imagine you are a detective trying to solve a mystery, but you have two mysterious machines, Machine A and Machine V.

• Machine A is a "Forward Machine." You put a piece of paper (data) in, and it spits out a result. You can see the result, but you don't know how the machine works inside.
• Machine V is a "Backward Machine." You put a result in, and it spits out a piece of paper. Again, you can see the output, but the inside is a black box.

In the world of medical imaging (like CT scans), these machines represent how X-rays travel through the body (Forward) and how we try to reconstruct the image from those X-rays (Backward). Ideally, these two machines should be perfect mirror images of each other. If you send a signal through A and then immediately through V, you should get exactly what you started with.

The Problem: In reality, these machines are often built by different people or using different math shortcuts. They aren't perfect mirrors. There is a "mismatch" or a "glitch" between them. The scientists want to know: How big is this glitch?

Mathematically, they want to calculate the "Operator Norm" of the difference between these two machines. But here's the catch:

1. They can't open the machines to see the gears (no access to the internal matrix).
2. They can't store a massive list of every possible input and output (the memory is too small).
3. They need a way to measure the glitch that gets more accurate the longer they run it.

The Solution: The "Blindfolded Hiker" Algorithm

The authors propose a clever, randomized method to measure this glitch. Think of it like a blindfolded hiker trying to find the highest peak in a foggy mountain range.

1. The Starting Point: The hiker (the algorithm) picks a random spot on the mountain (a random input vector) and checks the height (the output).
2. The "Adjoint" Twist: Usually, to climb a mountain efficiently, you need to know the slope (the gradient). But since the machines are black boxes, the hiker can't see the slope directly.
   • However, the hiker has a special trick: they can ask Machine A for the height, and they can ask Machine V for the "reverse" height. By comparing these two, they can estimate the direction of the steepest climb without seeing the map.
3. The Random Search: The hiker doesn't just walk straight up. They take a random step in a direction perpendicular to where they are standing.
   • They calculate the "best step size" (how far to walk) to maximize the difference they are measuring.
   • If the step makes the "glitch" look bigger, they keep it. If not, they adjust.
4. The Magic of Randomness: Because the hiker is taking random steps in all directions over and over, they are guaranteed to eventually stumble upon the exact direction where the glitch is the biggest.

Why This is Special

Most standard methods for finding the "biggest glitch" (like the Power Method) require you to know the exact blueprint of the machine (the adjoint matrix). If you don't have the blueprint, or if the machine is too big to fit in your computer's memory, those methods fail.

This new method is like a survivalist who can navigate a forest without a map or a compass, just by feeling the wind and the terrain.

• Memory Efficient: It only needs to remember two small pieces of paper (vectors) at a time, rather than a whole library of data.
• Guaranteed to Work: The paper proves mathematically that if you keep running this random search, you will almost certainly find the true size of the mismatch. It won't get stuck in a small valley; it will find the highest peak.

Real-World Application: The CT Scan

The authors tested this on Radon Transforms, which are the math behind CT scanners.

• In a CT scanner, the "Forward" part is the X-ray beam passing through you.
• The "Backward" part is the computer trying to build the 3D image of your bones.
• Sometimes, the software used to build the image isn't the perfect mathematical opposite of the software that simulates the X-ray. This causes blurry images or errors.

Using their new algorithm, the scientists could plug in the "Forward" code and the "Backward" code (without seeing the math inside) and calculate exactly how much they mismatched. They found that for some standard medical software, the mismatch was tiny (perfect), but for others, it was significant.

The Takeaway

This paper gives us a new tool to measure the "distance" between two complex, hidden systems. It's a stochastic (random) search method that acts like a determined explorer, using only the inputs and outputs of black-box machines to find the maximum error between them.

It's like trying to find the loudest echo in a cave by shouting random noises and listening carefully, rather than needing to see the cave's walls. It's efficient, it's smart, and it works even when you can't see the whole picture.

Technical Summary

Here is a detailed technical summary of the paper "Computing adjoint mismatch of linear maps" by Bresch et al.

1. Problem Statement

The paper addresses the computational challenge of estimating the operator norm (spectral norm) of the difference between two linear maps, $A$ and $V$, denoted as $\|A - V\|$. This problem arises in specific scenarios with restrictive constraints:

• Black-box Access: The map $A$ is only available as a forward operator (computing $Av$), while the map $V$ is only available via its adjoint (computing $V^*u$). The full matrix representations of $A$ and $V$ are not available.
• Memory Constraints: The dimensions of the spaces ($d$ for input, $m$ for output) are large, making it impossible to store dense matrices or a large number of vectors. The algorithm must have a storage complexity of $O(\max\{m, d\})$.
• Application Context: This is particularly relevant in Computerized Tomography (CT), where forward projection (Radon transform) and back-projection are often discretized independently. These discretizations frequently result in "adjoint mismatch," where the discrete back-projection is not the true adjoint of the forward projection. Quantifying $\|A - V\|$ is crucial for analyzing the convergence and error bounds of iterative reconstruction methods (e.g., Chambolle-Pock) that rely on these operators.

2. Methodology

The authors propose a stochastic algorithm (Algorithm 1) that generalizes the Rayleigh quotient maximization to find the largest singular value of the operator $L = A - V$.

Core Concept

The operator norm is defined as:
$$\|A - V\| = \max_{\|u\|=1, \|v\|=1} \langle u, (A - V)v \rangle$$
Since $V$ is only accessible via its adjoint, the term $\langle u, (A - V)v \rangle$ is computed as $\langle u, Av \rangle - \langle V^*u, v \rangle$.

Algorithmic Steps

1. Initialization: Randomly initialize unit vectors $u_0 \in \mathbb{R}^m$ and $v_0 \in \mathbb{R}^d$. Ensure the initial objective value is non-negative.
2. Random Search Directions: At each iteration $k$, sample random search directions $x_k$ and $w_k$ uniformly from the tangent spaces of the current vectors $v_k$ and $u_k$ (i.e., orthogonal to $v_k$ and $u_k$).
3. Optimal Step Size Calculation: Instead of using a fixed or heuristic step size, the algorithm analytically solves for the optimal step sizes $(\tau_k, \xi_k)$ that maximize the objective function along the search directions.
   • The objective function is reduced to a rational function of two variables.
   • The authors derive a closed-form solution for the maximizer of the squared objective, involving the calculation of specific inner products ($a_k, b_k, c_k, d_k$).
4. Update: The vectors are updated via:
   $$u_{k+1} = \frac{u_k + \tau_k w_k}{\|u_k + \tau_k w_k\|}, \quad v_{k+1} = \frac{v_k + \xi_k x_k}{\|v_k + \xi_k x_k\|}$$
5. Estimation: The estimate of the norm is updated as $a_k = \langle u_k, (A - V)v_k \rangle$.

Theoretical Properties

• Monotonicity: The sequence of objective values is strictly increasing (almost surely).
• Convergence: The sequence converges almost surely to the largest singular value $\sigma_1 = \|A - V\|$.
• Singular Vectors: The sequence of vector pairs $(u_k, v_k)$ converges almost surely to the corresponding left and right singular vectors associated with the largest singular value.
• Convergence Rate: The paper establishes a convergence rate of $O(1/n)$ for the probability of the error in the eigenvector equation exceeding a threshold. Numerical experiments suggest the practical convergence is often faster (nearly exponential).

3. Key Contributions

1. Novel Algorithm: The first method to compute $\|A - V\|$ under the specific "black-box" constraint where $A$ is forward-only and $V$ is adjoint-only, without requiring matrix assembly.
2. Optimal Step Sizes: Unlike standard stochastic gradient methods that use diminishing or fixed step sizes, this method computes exact optimal step sizes for the local search direction at every iteration, significantly accelerating convergence.
3. Minimal Storage: The algorithm requires storing only four vectors (two of dimension $d$, two of dimension $m$), satisfying the $O(\max\{m, d\})$ storage constraint, which is optimal for this problem class.
4. Rigorous Convergence Analysis: The authors provide a comprehensive proof of almost sure convergence to the operator norm and the corresponding singular vectors, handling cases where singular values have multiplicities.
5. Stopping Criteria: The paper proposes a practical stopping criterion based on the convergence of the parameters $b_k$ and $c_k$ to zero.

4. Results and Numerical Experiments

The authors validate the method through several experiments:

• Synthetic Data: Tests on random Gaussian matrices of varying dimensions ($m \times d$) show that the algorithm converges to the true norm. The relative error decreases as iterations increase.
• Comparison with [4]: The method is compared to a previous algorithm (from reference [4]) designed for computing $\|A\|$ (where $V=0$).
  • Finding: The proposed algorithm (Algorithm 1) is slightly slower in convergence rate than the specialized method in [4] when $V=0$, but it is the only method capable of handling the general case where $V \neq 0$ and $V^*$ is the only available form of $V$.
• Radon Transform (Astra Toolbox): The method was applied to check the adjointness of forward and back-projection operators in the Astra toolbox (a standard CT software).
  • Result: The algorithm confirmed that specific implementations (Line model, Ray model, Joseph method) are effectively adjoint (norm difference $\approx 10^{-9}$).
  • Contrast: When applied to standard MATLAB radon and iradon functions, the algorithm detected a significant mismatch (relative norm difference $> 0.1$), demonstrating its utility in identifying discretization errors in real-world tomographic systems.

5. Significance

This work is significant for the fields of inverse problems, medical imaging, and numerical linear algebra:

• Practical Utility: It enables the quantification of "adjoint mismatch" in large-scale imaging problems where storing the full system matrix is impossible. This is critical for ensuring the stability and accuracy of iterative reconstruction algorithms.
• Theoretical Advancement: It extends the theory of stochastic power methods to a "mixed" operator setting (one forward, one adjoint) and provides a rigorous convergence analysis for a method using optimal step sizes in a stochastic setting.
• Black-Box Flexibility: The approach allows researchers to verify the mathematical consistency of complex, proprietary, or highly optimized black-box operators (common in industrial CT software) without needing access to their internal matrix structures.

In summary, the paper provides a robust, memory-efficient, and theoretically sound tool for diagnosing and quantifying the discrepancy between forward and adjoint operators in high-dimensional linear systems.