The normal-inverse-Wishart (NIW) distribution is commonly used as a prior distribution for the mean and covariance parameters of a multivariate normal distribution. The family of NIW distributions is also a minimal exponential family. In this short note we describe a convergent procedure for converting from mean parameters to natural parameters in the NIW family, or -- equivalently -- for performing maximum likelihood estimation of the natural parameters given observed sufficient statistics. This is needed, for example, when using a NIW base family in expectation propagation
The Lippmann--Schwinger--Lanczos (LSL) algorithm has recently been shown to provide an efficient tool for imaging and direct inversion of synthetic aperture radar data in multi-scattering environments [17], where the data set is limited to the monostatic, a.k.a. single input/single output (SISO) measurements. The approach is based on constructing data-driven estimates of internal fields via a reduced-order model (ROM) framework and then plugging them into the Lippmann-Schwinger integral equation. However, the approximations of the internal solutions may have more error due to missing the off diagonal elements of the multiple input/multiple output (MIMO) matrix valued transfer function. This, in turn, may result in multiple echoes in the image. Here we present a ROM-based data completion algorithm to mitigate this problem. First, we apply the LSL algorithm to the SISO data as in [17] to obtain approximate reconstructions as well as the estimate of internal field. Next, we use these estimates to calculate a forward Lippmann-Schwinger integral to populate the missing off-diagonal data (the lifting step). Finally, to update the reconstructions, we solve the Lippmann-Schwinger equation using the original SISO data, where the internal fields are constructed from the lifted MIMO data. The steps of obtaining the approximate reconstructions and internal fields and populating the missing MIMO data entries can be repeated for complex models to improve the images even further. Efficiency of the proposed approach is demonstrated on 2D and 2.5D numerical examples, where we see reconstructions are improved substantially.
Predicting quantum operator matrices such as Hamiltonian, overlap, and density matrices in the density functional theory (DFT) framework is crucial for understanding material properties. Current methods often focus on individual operators and struggle with efficiency and scalability for large systems. Here we introduce a novel deep learning model, SLEM (Strictly Localized Equivariant Message-passing) for predicting multiple quantum operators, that achieves state-of-the-art accuracy while dramatically improving computational efficiency. SLEM's key innovation is its strict locality-based design, constructing local, equivariant representations for quantum tensors while preserving physical symmetries. This enables complex many-body dependence without expanding the effective receptive field, leading to superior data efficiency and transferability. Using an innovative SO(2) convolution technique, SLEM reduces the computational complexity of high-order tensor products and is therefore capable of handling systems requiring the $f$ and $g$ orbitals in their basis sets. We demonstrate SLEM's capabilities across diverse 2D and 3D materials, achieving high accuracy even with limited training data. SLEM's design facilitates efficient parallelization, potentially extending DFT simulations to systems with device-level sizes, opening new possibilities for large-scale quantum simulations and high-throughput materials discovery.
We investigate the set of invariant idempotent probabilities for countable idempotent iterated function systems (IFS) defined in compact metric spaces. We demonstrate that, with constant weights, there exists a unique invariant idempotent probability. Utilizing Secelean's approach to countable IFSs, we introduce partially finite idempotent IFSs and prove that the sequence of invariant idempotent measures for these systems converges to the invariant measure of the original countable IFS. We then apply these results to approximate such measures with discrete systems, producing, in the one-dimensional case, data series whose Higuchi fractal dimension can be calculated. Finally, we provide numerical approximations for two-dimensional cases and discuss the application of generalized Higuchi dimensions in these scenarios.
Finding an approximation of the inverse of the covariance matrix, also known as precision matrix, of a random vector with empirical data is widely discussed in finance and engineering. In data-driven problems, empirical data may be ``contaminated''. This raises the question as to whether the approximate precision matrix is reliable from a statistical point of view. In this paper, we concentrate on a much-noticed sparse estimator of the precision matrix and investigate the issue from the perspective of distributional stability. Specifically, we derive an explicit local Lipschitz bound for the distance between the distributions of the sparse estimator under two different distributions (regarded as the true data distribution and the distribution of ``contaminated'' data). The distance is measured by the Kantorovich metric on the set of all probability measures on a matrix space. We also present analogous results for the standard estimators of the covariance matrix and its eigenvalues. Furthermore, we discuss two applications and conduct some numerical experiments.
We propose an extremely versatile approach to address a large family of matrix nearness problems, possibly with additional linear constraints. Our method is based on splitting a matrix nearness problem into two nested optimization problems, of which the inner one can be solved either exactly or cheaply, while the outer one can be recast as an unconstrained optimization task over a smooth real Riemannian manifold. We observe that this paradigm applies to many matrix nearness problems of practical interest appearing in the literature, thus revealing that they are equivalent in this sense to a Riemannian optimization problem. We also show that the objective function to be minimized on the Riemannian manifold can be discontinuous, thus requiring regularization techniques, and we give conditions for this to happen. Finally, we demonstrate the practical applicability of our method by implementing it for a number of matrix nearness problems that are relevant for applications and are currently considered very demanding in practice. Extensive numerical experiments demonstrate that our method often greatly outperforms its predecessors, including algorithms specifically designed for those particular problems.
Anderson acceleration (AA) is widely used for accelerating the convergence of an underlying fixed-point iteration $\bm{x}_{k+1} = \bm{q}( \bm{x}_{k} )$, $k = 0, 1, \ldots$, with $\bm{x}_k \in \mathbb{R}^n$, $\bm{q} \colon \mathbb{R}^n \to \mathbb{R}^n$. Despite AA's widespread use, relatively little is understood theoretically about the extent to which it may accelerate the underlying fixed-point iteration. To this end, we analyze a restarted variant of AA with a restart size of one, a method closely related to GMRES(1). We consider the case of $\bm{q}( \bm{x} ) = M \bm{x} + \bm{b}$ with matrix $M \in \mathbb{R}^{n \times n}$ either symmetric or skew-symmetric. For both classes of $M$ we compute the worst-case root-average asymptotic convergence factor of the AA method, partially relying on conjecture in the symmetric setting, proving that it is strictly smaller than that of the underlying fixed-point iteration. For symmetric $M$, we show that the AA residual iteration corresponds to a fixed-point iteration for solving an eigenvector-dependent nonlinear eigenvalue problem (NEPv), and we show how this can result in the convergence factor strongly depending on the initial iterate, which we quantify exactly in certain special cases. Conversely, for skew-symmetric $M$ we show that the AA residual iteration is closely related to a power iteration for $M$, and how this results in the convergence factor being independent of the initial iterate. Supporting numerical results are given, which also indicate the theory is applicable to the more general setting of nonlinear $\bm{q}$ with Jacobian at the fixed point that is symmetric or skew symmetric.
An equivalence relation can be constructed from a given (homogeneous, binary) relation in two steps: first, construct the smallest reflexive and transitive relation containing the given relation (the "star" of the relation) and, second, construct the largest symmetric relation that is included in the result of the first step. The fact that the final result is also reflexive and transitive (as well as symmetric), and thus an equivalence relation, is not immediately obvious, although straightforward to prove. Rather than prove that the defining properties of reflexivity and transitivity are satisfied, we establish reflexivity and transitivity \emph{constructively} by exhibiting a particular starth root -- in a way that emphasises the creative process in its construction. The constructed starth root is fundamental to algorithms that determinethe strongly connected components of a graph as well as the decomposition of a graph into its strongly connected components together with an acyclic graph connecting such components.
In this paper, we effectively solve the inverse source problem of the fractional Poisson equation using MC-fPINNs. We construct two neural networks $ u_{NN}(x;\theta )$ and $f_{NN}(x;\psi)$ to approximate the solution $u^{*}(x)$ and the forcing term $f^{*}(x)$ of the fractional Poisson equation. To optimize these two neural networks, we use the Monte Carlo sampling method mentioned in MC-fPINNs and define a new loss function combining measurement data and the underlying physical model. Meanwhile, we present a comprehensive error analysis for this method, along with a prior rule to select the appropriate parameters of neural networks. Several numerical examples are given to demonstrate the great precision and robustness of this method in solving high-dimensional problems up to 10D, with various fractional order $\alpha$ and different noise levels of the measurement data ranging from 1$\%$ to 10$\%$.
This paper addresses second-order stochastic optimization for estimating the minimizer of a convex function written as an expectation. A direct recursive estimation technique for the inverse Hessian matrix using a Robbins-Monro procedure is introduced. This approach enables to drastically reduces computational complexity. Above all, it allows to develop universal stochastic Newton methods and investigate the asymptotic efficiency of the proposed approach. This work so expands the application scope of secondorder algorithms in stochastic optimization.
A new, more efficient, numerical method for the SDOF problem is presented. Its construction is based on the weak form of the equation of motion, as obtained in part I of the paper, using piece-wise polynomial functions as interpolation functions. The approximation rate can be arbitrarily high, proportional to the degree of the interpolation functions, tempered only by numerical instability. Moreover, the mechanical energy of the system is conserved. Consequently, all significant drawbacks of existing algorithms, such as the limitations imposed by the Dahlqvist Barrier theorem and the need for introduction of numerical damping, have been overcome.