FDGI 0.9

A package for frequency-domain ghost imaging

Spooktroscopy: Frequency-Domain Ghost Imaging

Target Problem

Solve in the least-square way, under regularizations. A, G are matrices acting on the two indices of X, i.e.

G is optional, by default (G=None) it is identity, in which case this is the conventional Spooktroscopy, i.e. to solve under regularizations.

The reason for G is to accommodate the combination with pBASEX, in which case this is solving the two linear inversions in one step.

Key Advantages

The key advantages of this package are

1. Efficient optimization: Contraction over shots is decoupled from optimization. It is recommended to input precontracted results when instantiating, or to save the caches using method save_prectr of a solver created with raw inputs. Once instantiated, the precontracted results are cached to be reused every time solving with a different hyperparameter. See spook.contraction_utils.adaptive_contraction .
2. Support dimension reduction on the dependent variable B, through basis functions in G. In that case, it is also recommended to contract (B,G) over the dependent variable space (index q) prior to instantiating a spook solver.
3. Support multiple combinations of regularizations. See Solvers .
4. Support time-dependent measurement (Beta): when each entry in the raw input A is a pair of (photon spectrum, delay bin), index w is the flattened axis of (, ). In this case, the third smoothness hyperparameter is for the delay axis.

At the very bottom level, this package depends on either OSQP to solve a quadratic programming or LAPACK gesv through numpy.linalg.solve .

Regularizations

Common regularizations are the following three types, all of which optional, depending on what a prior knowledge one wants to enforce on the problem solving.

1. Nonnegativity: To constrain
2. Sparsity: To penalize on or
3. Smoothness: To penalize on roughness of X , along the two indices, independently. For -axis, which is the photon energy axis, the form is fixed where is the laplacian. Roughness along the second axis of X is customizable through parameter Bsmoother, which by default is laplacian squared too.

Sparsity and Smoothness are enforced through penalties in the total obejctive function, and the penalties are weighted by hyperparameters lsparse and lsmooth. lsmooth is a 2-tuple that weight roughness penalty along the two axes of X respectively. The hyperparameters can be passed in during instantiation and also updated afterwards. It is recommended to call method getXopt with the hyperparameter(s) to be updated, because it will update, solve, and return the optimal X in one step. Calling solve with the hyperparameter(s) to be updated and then calling getXopt() without input is effectively the same, and the problem will be solved once as long as there is no update.

Solvers

Different combinations of regularizations can lead to different forms of objective function. Solvers in package always formalize the specific problem into either a Quadratic Programming or a linear equation. Examples can be found in unit tests

Nonnegativity	Sparsity	Smoothness	Solver	Notes
True	L1 or False	Quadratic	SpookPosL1	This solver can serve tasks like in Li et al
True	L2 squared	Quadratic	SpookPosL2
False	L2 squared or False	Quadratic	SpookLinSolve	This solver is so far the work-horse for SpookVMI
False	L1	Quadratic	SpookL1

Quadratic Programming

For cases where it can be formalized into a Quadratic Programming , OSQP does the job. Thus the root numerical method is alternating direction method of multipliers (ADMM). Looking into the solver settings of OSQP is always encouraged, but the default settings usually work fine for spook . If one needs to pass in settings, the OSQP solver is SpookQPBase._prob .

Linear Equation

A rare case that it can be formalized into a linear equation is the third line in the table above: no nonnegativity constraint, and the sparsity is L2 norm squared. This is implemented in SpookLinSolve , which calls numpy.linalg.solve or scipy.sparse.linalg.spsolve .

Normalization Convention

The entries in are preferred to be on the order of unity, because regularization-related quadratic form matrices have their entries around unity. The scale factors are set as

where are the dimensions along w-axis and q-axis, respectively. is an accessible property of the solver AGscale. To normalize or not is controlled by parameter pre_normalize in instanciation.

By default pre_normalize=True, i.e. self._AtA =, self._GtG =, and self._Bcontracted = . In this case, the direct solution self.res is scaled as , but in getXopt the final result of is scaled back before returned.

The entries in B are not always accessible, because of the option to pass in precontracted results and in mode='contracted'. Therefore B is not normalized.

Unit Tests

unittest/testPosL1.py is a good example to play with SpookPosL1. unittest/testL1L2.ipynb include good examples to play with SpookPosL1, SpookPosL2, and SpookL1.

Dependencies

See requirements.txt

Acknowledgement

This work was supported by the U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences (BES), Chemical Sciences, Geosciences, and Biosciences Division (CSGB).