Fast Fusion of Multi-Band Images Based on Solving a Sylvester Equation

Qi Wei; Nicolas Dobigeon; Jean-Yves Tourneret

Source-linked AI summary

Fast Fusion of Multi-Band Images Based on Solving a Sylvester Equation

Qi Wei, Nicolas Dobigeon, Jean-Yves Tourneret

arXiv:1502.03121v1Feb 10, 2015cs.CV

TL;DR

Multi-band fusion must reconstruct a high-spatial and high-spectral image from complementary observations, while existing estimation methods can be computationally expensive. The paper derives an explicit Sylvester-equation solution using circulant blurring and downsampling structure, and reports competitive performance with substantially lower computational complexity. It also extends the approach to Bayesian estimators using ADMM and block coordinate descent.

Problem
Reconstructing a high-spatial and high-spectral image from complementary observations is computationally challenging, especially for large images and downsampling models.
Method
The method solves the fusion problem through an explicit Sylvester-equation solution and combines it with ADMM and block coordinate descent for Bayesian estimation.
Results
The proposed method compares competitively with ADMM-based methods while significantly reducing computational complexity.
Takeaways & Limitations
The closed-form Sylvester solution makes fast fusion applicable to multispectral–hyperspectral and panchromatic–hyperspectral image fusion.
Takeaways & Limitations
Only the Gaussian prior of is considered in the reported experiment because of space limitations.

Abstract

from arXiv · show

This paper proposes a fast multi-band image fusion algorithm, which combines a high-spatial low-spectral resolution image and a low-spatial high-spectral resolution image. The well admitted forward model is explored to form the likelihoods of the observations. Maximizing the likelihoods leads to solving a Sylvester equation. By exploiting the properties of the circulant and downsampling matrices associated with the fusion problem, a closed-form solution for the corresponding Sylvester equation is obtained explicitly, getting rid of any iterative update step. Coupled with the alternating direction method of multipliers and the block coordinate descent method, the proposed algorithm can be easily generalized to incorporate prior information for the fusion problem, allowing a Bayesian estimator. Simulation results show that the proposed algorithm achieves the same performance as existing algorithms with the advantage of significantly decreasing the computational complexity of these algorithms.

I. INTRODUCTION

Multi-band fusion reconstructs a high-spatial and high-spectral image from complementary degraded observations, but large-scale estimation is computationally challenging. The paper addresses this challenge with an explicit Sylvester-equation solution that supports efficient ML and Bayesian fusion.

Motivation: Multi-band images trade spatial resolution against spectral sensitivity, motivating fusion of complementary observations.Hyperspectral images provide hundreds of bands but lower spatial resolution than multispectral or panchromatic images.
Problem formulation: The fusion task estimates an unknown high-spatial and high-spectral image from two degraded data cubes.The observations undergo complementary spatial and spectral transformations.
Computational challenge: Large image dimensions and downsampling make ML and Bayesian estimation difficult by preventing direct Fourier diagonalization of the blurring operator.Existing computational strategies include approximate estimators and computationally expensive MCMC methods.
Proposed approach: The proposed method explicitly solves the Sylvester equation underlying the fusion problem and computes the ML estimator efficiently.The approach avoids the iterative methods used by preceding approaches.
Bayesian extension: The method generalizes to Bayesian fusion when appropriate prior information is incorporated.The paper describes this generalization as computationally efficient within a Bayesian framework.

C. Paper organization

The paper formulates fusion under a linear degradation model with sensor-dependent Gaussian noise and separates optimization into ML and Bayesian estimation. A low-dimensional subspace representation reduces the spectral parameterization.

Paper organization: The paper studies ML fusion without regularization and later generalizes the method to Bayesian estimators with various priors.The paper organization assigns these topics to Sections II–IV and evaluates accuracy and numerical efficiency in Section V.
Likelihood model: The two observed images are conditionally independent given the unobserved scene because their sensor-dependent noise matrices are assumed statistically independent.This independence determines the joint likelihood formulation.
Subspace representation: The target image is represented as X = HU, where U contains its projection onto a lower-dimensional spectral subspace.The subspace dimension is smaller than the number of hyperspectral bands because adjacent bands are highly correlated.

U L(U)

Maximizing the likelihood produces a generalized Sylvester equation for the target subspace coefficients. Circulant blurring and structured downsampling enable an analytical solution under stated assumptions.

Analytical solution: The paper proves that minimizing the likelihood with respect to U can be solved analytically without iterative optimization or Monte Carlo methods.The resulting closed-form solution is presented as the basis of the fast algorithm.
Sylvester formulation: The ML optimization with respect to U becomes a generalized Sylvester matrix equation after manipulating the normal equations.A unique solution requires that arbitrary sums of eigenvalues from the two coefficient matrices are nonzero.
Blurring assumption: Assuming a circulant blurring matrix allows B to be decomposed through the discrete Fourier transform into a diagonal matrix D.The conjugate blurring matrix has the corresponding decomposition using D*.
Downsampling assumption: The decimation matrix represents uniform downsampling, while its conjugate transpose performs zero interpolation.The relation S^T S = I_m supports the structured treatment of sampling in the derivation.

B. Existence of a solution

The Sylvester equation has a unique solution when the eigenvalue sums of C1 and C2 are nonzero, but this requires C1 to be nonsingular. The paper exploits the fusion-specific matrix structure to derive an analytically computable closed-form solution instead of using a costly general-purpose solver.

C1 is diagonalizable with non-negative eigenvalues because it is formed from positive-definite matrix factors.
A unique solution of (7) follows when C1 is invertible, so every sum of eigenvalues of C1 and C2 is greater than zero.
C1 invertibility is not guaranteed by the dimensions and forms of H and L; for nλ < emλ, L LH is rank deficient.
When C1 is singular, regularization or prior information is needed to ensure a unique solution, while the section focuses on the nonsingular case.
The classical Bartels-Stewart approach is impractical because C2 = BH is n × n and its Schur decomposition costs O(n^3) operations.
The proposed strategy exploits the specific structures of C1 and C2 to obtain an analytic Sylvester-equation solution, computed block-by-block with trivial diagonal-matrix inversion.

12 Set X = HQ ¯UPD−1FH; Output: X

The dominant computational costs arise from FFT-based matrix operations, while frequency-domain calculations reduce other operations to near-linear or structured updates.

The most computationally expensive steps compute D and ¯C3 using FFT and iFFT operations.
Decomposing B = FDFH, equivalently transforming the blurring kernel, has complexity O(n log n).
The calculations involving L YLF require one FFT operation each, while other computations occur in the frequency domain.
Multiplication by DP−1 costs O(n) because D is diagonal and P−1 reduces to block shifting and addition.
The computation of C3BFP−1 has total complexity of order O(n · max {log n, em2.

IV. GENERALIZATION TO BAYESIAN ESTIMATORS

The Bayesian generalization introduces priors when the Sylvester equation may lack a unique solution, and combines closed-form Sylvester solving with Gaussian and non-Gaussian regularization strategies.

When the MS-band count is smaller than the subspace dimension, a prior or regularization is required to ensure a unique Sylvester-equation solution.
Combining the likelihoods with p(U) yields a posterior, whose mode is the MAP estimator.
The explicit Sylvester-equation solution can compute MAP estimators for the considered prior distributions.
Gaussian priors: For a matrix normal prior, covariance Σ models HS-band correlations and controls distance from the mean µ.
Non-Gaussian priors: Non-Gaussian priors split the objective into likelihood and regularization terms and use ADMM with a Gaussian-prior surrogate and proximity operator.
Non-Gaussian priors: The U update can use the analytical solver, with per-iteration complexity O(n log n) from FFT and iFFT steps.

2) Solution in frequency domain:

The frequency-domain formulation embeds the Sylvester solver within ADMM or block coordinate descent, extending the closed-form approach to hierarchical and non-Gaussian Bayesian models.

The optimization can be solved by alternating updates of U and hyperparameters Φ under a hyperprior p(Φ).
The method is applied to generalized Gaussian, ℓ1 or ℓ0, ℓ2 smooth, and TV regularizations investigated for multi-band fusion.
The U update uses the closed-form solver, while the Φ update depends on the hyperprior form.
SE-BCD denotes the resulting SE-within-block-coordinate-descent method.
The approach treats the non-Gaussian-prior strategy as a special case of hierarchical updating.

A. Fusion Quality Metrics

Fusion quality is evaluated with five image-quality measures across several prior models and experiments; the proposed methods generally match competing methods while reducing computation.

A. Fusion Quality Metrics: Five measures evaluate fusion quality: RSNR, SAM, UIQI, ERGAS, and DD.
A. Fusion Quality Metrics: Larger RSNR and UIQI indicate better fusion, whereas smaller SAM, ERGAS, and DD indicate better fusion.
B. Fusion of HS and MS images: In the naive Gaussian-prior experiment, estimated images are visually similar and quality measures confirm similar performance.
B. Fusion of HS and MS images: The proposed computational time is reduced by factors larger than 100 for supervised fusion and 50 for unsupervised fusion.
B. Fusion of HS and MS images: For sparse representation, SE-BCD has similar performance to ADMM-BCD but significantly reduces complexity.
B. Fusion of HS and MS images: For TV regularization, SE-ADMM produces similar fusion results, reduced computational time, faster initial convergence, and smoother convergence.

C. Hyperspectral Pansharpening

The hyperspectral pansharpening experiments apply the proposed Sylvester-equation-based fusion method to simulated and real hyperspectral data. The method produces similar fusion quality to ADMM-based alternatives while substantially reducing computational complexity and supporting Bayesian extensions.

Experimental setting: Hyperspectral pansharpening extends conventional pansharpening by fusing hyperspectral and panchromatic images for remote sensing.The experiments use an AVIRIS reference image from Moffett Field, reduced from 224 to 176 bands after removing water-vapor absorption bands.
Experimental setting: The simulated HS and PAN observations are generated from a reference image with averaging, realistic spectral responses, and additive Gaussian noise.The HS image has SNRs of 35 dB for the first 126 bands and 30 dB for the remaining bands.
Methods and evaluation: A Gaussian prior regularizes the ill-posed supervised pansharpening problem, and the Sylvester-equation method is compared with ADMM.The comparison uses fixed hyperparameters, with visual results in Fig. 3 and quantitative results in Table III.
Results: The proposed SE-based method provides similar qualitative and quantitative fusion results to ADMM while requiring significantly lower computational complexity.The comparison evaluates multiple fusion-quality measures and computational efficiency.
Method scope: The explicit Sylvester-equation solution applies to multispectral–hyperspectral and panchromatic–hyperspectral fusion and can be generalized to Bayesian estimators with priors.The Bayesian generalization uses alternating direction method of multipliers and block coordinate descent.
Future work: Future work includes learning the subspace transform matrix H and implementing the fusion scheme on real datasets.These directions define the stated scope boundary of the reported development.

APPENDIX A

Appendix A develops matrix properties and lemmas used to analyze the Sylvester-equation solution. It establishes positivity, diagonalizability, and a blockwise procedure for obtaining solution elements.

Matrix properties: Positive-definite and positive-semidefinite matrix factors establish invertibility and nonnegative eigenvalues for the matrices used in the derivation.The argument also uses similarity to transfer eigenvalue and diagonalizability properties.
Matrix properties: The appendix shows that the relevant product is diagonalizable by relating it to a symmetric or Hermitian matrix.Similar matrices share eigenvalues, supporting the spectral analysis.
Lemma 3: Lemma 3 supplies an equality involving the matrices defined in Section III-A, the all-ones matrix J_d, and the identity matrix I_m.The proof is deferred to Appendix D.
Solution construction: Once the first-column entries are known, the remaining entries are obtained for l = 1, · · ·, emλ and j = 2, · · ·, d.This completes the blockwise construction described in the appendix.

APPENDIX D

Appendix D analyzes the discrete Fourier and elementary matrices used in the fusion derivation. It defines their structure and records properties such as idempotence.

DFT representation: The n-dimensional DFT matrix F is written explicitly using powers of a primitive nth root of unity.The root is defined through ω and i = √−1.
Elementary matrices: The matrix E_i is defined as an n × n matrix with a single nonzero entry equal to 1 at row i and column i.This gives an elementary selector-like matrix for the derivation.
Elementary matrices: Each E_i is idempotent, satisfying E_i = E_i^2.The appendix uses this elementary-matrix property in its subsequent identities.
DFT representation: The appendix specifies rows of F and the zero vector, then derives related matrix expressions by direct computation.These expressions support the subsequent structured matrix manipulations.
Block structure: The displayed structured matrices contain repeated zero blocks and identity blocks in the decomposition used by the proof.The block layout is part of the appendix’s explicit matrix construction.