In this paper, we show how informativity and identifiability for networks of dynamical systems can be investigated using Gröbner bases. We provide a sufficient condition for informativity in terms of positive definiteness of the spectrum of external signals and full generic rank of the transfer function relating the external signals to the inputs of the predictor. Moreover, we show how generic local network identifiability can be investigated by computing the dimension of the fiber associated with the closed loop transfer function from external measurable signals to the measured outputs.

In practical applications, measured time series are always associated with a certain amount of noise. This causes the resulting numerical derivatives to be associated with significant high-frequency noise and integrals with low-frequency noise, requiring special treatment of the numerical procedures to compute these quantities. To address this challenge, this paper presents a Gaussian process (GP) regression approach for the differentiation and integration of one-dimensional signals such as time series. Modeling the underlying time series as a GP can mitigate noise amplification and provides a quantitative measure of the uncertainty of the computed derivative or integral. The state space form of GPs, combined with the Kalman filter and Rauch-Tung-Striebel smoother, is used herein to achieve computational efficiency in the regression process. Two numerical examples are provided including a Duffing oscillator and a Lorenz system, to showcase the performance of the presented method in a simplified scenario. A numerical application consisting of a 2D offshore wind turbine subjected to wave loads is then shown, in which the noisy acceleration response measurements are used to estimate the displacement response. This scenario is particularly relevant in the context of structural health monitoring (SHM), in which displacement responses can be very informative on the structural performance and internal stresses, but they are difficult and expensive to measure. The performance of the proposed methodology is then compared with other standard methods. The illustrative examples make use of the Mat & eacute;rn kernel, while the formulation for other popular kernels is also presented.

In this work, we study the well-posedness of certain sparse regularized linear regression problems, i.e., the existence, uniqueness and continuity of the solution map with respect to the data. We focus on regularization functions that are convex piecewise linear, i.e., whose epigraph is polyhedral. This includes total variation on graphs and polyhedral constraints. We provide a geometric framework for these functions based on their connection to polyhedral sets and apply this to the study of the well-posedness of the corresponding sparse regularized linear regression problem. Particularly, we provide geometric conditions for well-posedness of the regression problem, compare these conditions to those for smooth regularization, and show the computational difficulty of verifying these conditions.

This article proposes a scalable, matrix-free approach to kernel-based regularization for finite impulse response estimation. Our methodology is based on Bayesian optimization, a gradient-free global optimization strategy that can handle noisy objective function evaluations and is suitable for indirect or stochastic hyperparameter estimation. A key assumption is that the kernel function only has a small number of parameters and that matrix-vector products with the kernel matrix are inexpensive. Popular kernels, such as the tuned-correlated (TC), diagonal-correlated (DC), and stable spline (SS) kernels, meet these criteria. To achieve scalability, we simultaneously exploit the structure of both the kernel matrix and the regressor matrix, which reduces memory requirements and computational cost. Our empirical results, which are based on randomly generated systems, demonstrate that the computational cost scales approximately linearly with the model order for the TC, DC, and SS kernels.

This paper explores the problem of finding a rank-structured approximation to a given Hermitian positive definite matrix. We first study a fundamental matrix nearness problem which seeks an approximation whose off-diagonal block must have low rank. Our investigation focuses on three objective functions: the spectral norm, Frobenius divergence, and log-determinant divergence. We derive closed-form solutions for the fundamental nearness problem, extending our analysis to include problems with a range restrictions on the off-diagonal block. Building upon these results, we propose a practical, scalable algorithm to compute a positive definite hierarchical off-diagonal low-rank approximation to a Hermitian positive definite matrix. Through numerical experiments, we explore the use of our approximation as a general-purpose preconditioner for the conjugate gradient method. Our results demonstrate that the log-determinant divergence-based approximation exhibits good performance across a diverse set of test problems.

In this experimental work, we present a general framework based on the Bregman log determinant divergence for preconditioning Hermitian positive definite linear systems. We explore this divergence as a measure of discrepancy between a preconditioner and a matrix. Given an approximate factorization of a given matrix, the proposed framework informs the construction of a low-rank approximation of the typically indefinite factorization error. The resulting preconditioner is therefore a sum of a Hermitian positive definite matrix given by an approximate factorization plus a low-rank matrix. Notably, the low-rank term is not generally obtained as a truncated singular value decomposition (TSVD). This framework leads to a new truncation where principal directions are not based on the magnitude of the singular values, and we prove that such truncations are minimizers of the aforementioned divergence. We present several numerical examples showing that the proposed preconditioner can reduce the number of PCG iterations compared to a preconditioner constructed using a TSVD for the same rank. We also propose a heuristic to approximate the proposed preconditioner in the case where exact truncations cannot be computed explicitly (e.g., in a large-scale setting) and demonstrate its effectiveness over TSVD-based approaches.

Maximum likelihood estimation is effective for identifying dynamical systems, but applying it to large networks becomes computationally prohibitive. This paper introduces a maximum likelihood estimation method that enables identification of sub-networks within complex interconnected systems without estimating the entire network. The key insight is that under specific topological conditions, a sub-network's parameters can be estimated using only local measurements: signals within the target sub-network and those in the directly connected to the so-called separator sub-network. This approach significantly reduces computational complexity while enhancing privacy by eliminating the need to share sensitive internal data across organizational boundaries. We establish theoretical conditions for network separability, derive the probability density function for the sub-network, and demonstrate the method's effectiveness through numerical examples.

Numerically efficient and stable implementation of algorithms is essential for the kernel-based regularized system identification in practice. The state of art algorithms explore the semiseparable structure of the kernel and are based on the generator representation of the kernel matrix. However, as will be shown from both the theory and the practice, the algorithms based on the generator representation are sometimes numerically unstable, and thus limits its application in practice. In this paper, we aim to address this issue, and we consider the alternative Givens-vector representation of semiseparable kernels instead, which is numerically more stable but often much harder to derive. In particular, we derive the Givens-vector representation of some widely used kernel matrices. Then, we design algorithms based on the Givens-vector representation. Monte Carlo simulations show that the proposed algorithms admit the same order of computational complexity as the state-of-the-art ones based on generator representation, but with more stable and accurate implementation.

This paper investigates maximum likelihood estimation for direct system identification in networks of dynamical systems. We establish that the proposed approach is both consistent and efficient. In addition, it is more generally applicable than existing methods, since it can be employed even when measurements are unavailable for all network nodes, provided that network identifiability is satisfied. Finally, we demonstrate that the maximum likelihood problem can be formulated without relying on a predictor, which is key to achieving computationally efficient numerical solutions.

Numerically efficient and stable algorithms are essential for kernel-based regularized system identification. The state of art algorithms exploit the semiseparable structure of the kernel and are based on the generator representation of the kernel matrix. However, as will be shown from both the theory and the practice, the algorithms based on the generator representation are sometimes numerically unstable, which limits their application in practice. This paper aims to address this issue by deriving and exploiting an alternative Givens-vector representation of some widely used kernel matrices. Based on the Givens-vector representation, we derive algorithms that yield more accurate results than existing algorithms without sacrificing efficiency. We demonstrate their usage for the kernel-based regularized system identification. Monte Carlo simulations show that the proposed algorithms admit the same order of computational complexity as the state-of-the-art ones based on generator representation, but without issues with numerical stability.

This paper studies two classes of sampling methods for the solution of inverse problems, namely Randomize-Then-Optimize (RTO), which is rooted in sensitivity analysis, and Langevin methods, which are rooted in the Bayesian framework. The two classes of methods correspond to different assumptions and yield samples from different target distributions. We highlight the main conceptual and theoretical differences between the two approaches and compare them from a practical point of view by tackling two classical inverse problems in imaging: deblurring and inpainting. We show that the choice of the sampling method has a significant impact on the quality of the reconstruction and that the RTO method is more robust to the choice of the parameters.

This paper introduces a novel direct approach to system identification of dynamic networks with missing data based on maximum likelihood estimation. Dynamic networks generally present a singular probability density function, which poses a challenge in the estimation of their parameters. By leveraging knowledge about the network's interconnections, we show that it is possible to transform the problem into a more tractable form by applying linear transformations. This results in a nonsingular probability density function, enabling the application of maximum likelihood estimation techniques. Our preliminary numerical results suggest that when combined with global optimization algorithms or a suitable initialization strategy, we are able to obtain a good estimate of the dynamics of the internal systems.

This paper presents a fast direct solver for the Combined Field Integral Equation using higher-order discretizations. By adopting higher-order polynomials with the Method of Moments, the number of unknowns is significantly reduced. The fast direct solver leverages the efficiency of the Multi Level Fast Multipole Method by combining it with randomized linear algebra to construct low-rank approximations in a H-2 format. The proposed method is fully error controllable and achieves a setup time with computational complexity of O(r(3) logN). Numerical results for the scattering problem of a sphere demonstrate high accuracy, and the efficiency is demonstrated on the NASA Almond.

We study a preconditioner for a Hermitian positive definite linear system, which is obtained as the solution of a matrix nearness problem based on the Bregman log determinant divergence. The preconditioner is of the form of a Hermitian positive definite matrix plus a lowrank matrix. For this choice of structure, the generalized eigenvalues of the preconditioned matrix are easily calculated, and we show under which conditions the preconditioner minimizes the \ell 2 condition number of the preconditioned matrix. We develop practical numerical approximations of the preconditioner based on the randomized singular value decomposition (SVD) and the Nystro"\m approximation and provide corresponding approximation results. Furthermore, we prove that the Nystro"\m approximation is in fact also a matrix approximation in a range-restricted Bregman divergence and establish several connections between this divergence and matrix nearness problems in different measures. Numerical examples are provided to support the theoretical results.

We introduce AdaSub, a stochastic optimization algorithm that computes a search direction based on second-order information in a low-dimensional subspace that is defined adaptively based on available current and past information. Compared to first-order methods, second-order methods exhibit better convergence characteristics, but the need to compute the Hessian matrix at each iteration results in excessive computational expenses, making them impractical. To address this issue, our approach enables the management of computational expenses and algorithm efficiency by enabling the selection of the subspace dimension for the search. Our code is freely available on GitHub, and our preliminary numerical results demonstrate that AdaSub surpasses popular stochastic optimizers in terms of time and number of iterations required to reach a given accuracy.

Spectral computed tomography has received considerable interest in recent years since spectral measurements contain much richer information about the object of interest. In spectral computed tomography, we are interested in the energy channel-wise reconstructions of the object. However, such reconstructions suffer from low signal-to-noise ratio and share the challenges of conventional low-dose computed tomography such as ring artifacts. Ring artifacts arise from errors in the flat-field correction and can significantly degrade the quality of the reconstruction. We propose an extended flat-field model that exploits high correlation in the spectral flat-fields to reduce ring artifacts in the channel-wise reconstructions. The extended model relies on the assumption that the spectral flat-fields can be well-approximated by a low-rank matrix. Our proposed model works directly on the spectral flat-fields and can be combined with any existing reconstruction model, e.g., filtered back projection and iterative methods. The proposed model is validated on a neutron data set. The results show that our method successfully diminishes ring artifacts and improves the quality of the reconstructions. Moreover, the results indicate that our method is robust; it only needs a single spectral flat-field image, whereas existing methods need multiple spectral flat-field images to reach a similar level of ring reduction.

Computed tomography is a method for synthesizing volumetric or cross-sectional images of an object from a collection of projections. Popular reconstruction methods for computed tomography are based on idealized models and assumptions that may not be valid in practice. One such assumption is that the exact projection geometry is known. The projection geometry describes the relative location of the radiation source, object, and detector for each projection. However, in practice, the geometric parameters used to describe the position and orientation of the radiation source, object, and detector are estimated quantities with uncertainty. A failure to accurately estimate the geometry may lead to reconstructions with severe misalignment artifacts that significantly decrease their scientific or diagnostic value. We propose a novel reconstruction method that jointly estimates the reconstruction and the projection geometry. The reconstruction method is based on a Bayesian approach that yields a point estimate for the reconstruction and geometric parameters and, in addition, provides valuable information regarding their uncertainty. This is achieved by approximately sampling from the joint posterior distribution of the reconstruction and projection geometry using a hierarchical Gibbs sampler. Using real tomographic data, we demonstrate that the proposed reconstruction method significantly reduces misalignment artifacts. Compared with two commonly used alignment methods, our proposed method achieves comparable or better results under challenging conditions.

We present an algorithm for computing the minimum-rank positive semidefinite completion of a sparse matrix with a chordal sparsity pattern.This problem is tractable, in contrast to the minimumrank positive semidefinite completion problem for general sparsity patterns.We also present a similar algorithm for the Euclidean distance matrix completion with minimum embedding dimension.The two algorithms use efficient recursions over a clique tree associated with the chordal sparsity pattern.As an application, we use the minimum-rank completion method as a rounding technique to convert the solution of a sparse semidefinite optimization problem with non-unique solutions to an optimal solution of lower rank.In experiments with semidefinite relaxations of optimal power flow problems, the minimum-rank completion often results in solutions of lower rank than the solutions computed by interior-point solvers.

Constraints are a natural choice for prior information in Bayesian inference. In various applications, the parameters of interest lie on the boundary of the constraint set. In this paper, we use a method that implicitly defines a constrained prior such that the posterior assigns positive probability to the boundary of the constraint set. We show that by projecting posterior mass onto the constraint set, we obtain a new posterior with a rich probabilistic structure on the boundary of that set. If the original posterior is a Gaussian, then such a projection can be done efficiently. We apply the method to Bayesian linear inverse problems, in which case samples can be obtained by repeatedly solving constrained least squares problems, similar to a MAP estimate, but with perturbations in the data. When combined into a Bayesian hierarchical model and the constraint set is a polyhedral cone, we can derive a Gibbs sampler to efficiently sample from the hierarchical model. To show the effect of projecting the posterior, we applied the method to deblurring and computed tomography examples.

This paper proposes a methodology for scalable kernel-based regularized system identification based on indirect methods. It leverages stochastic trace estimation methods and an iterative solver such as LSQR for the efficient evaluation of hyperparameter selection criteria. It also uses a derivative-free optimization approach to hyperparameter estimation, which avoids the need for computing gradients or Hessians of the objective function. Moreover, the method is matrix-free, which means it only relies on a matrix-vector oracle and exploits fast routines for various structured matrix-vector products. Our preliminary numerical experiments indicate that the methodologygy scales significantly better than direct methods, especially when dealing with large datasets and slowly decaying impulse responses.