Existence and Uniqueness of Perturbation Solutions to DSGE Models

We prove that standard regularity and saddle stability assumptions for linear approximations are sufficient to guarantee the existence of a unique solution for all undetermined coefficients of nonlinear perturbations of arbitrary order to discrete time DSGE models. We derive the perturbation using a matrix calculus that preserves linear algebraic structures to arbitrary orders of derivatives, enabling the direct application of theorems from matrix analysis to prove our main result. As a consequence, we provide insight into several invertibility assumptions from linear solution methods, prove that the local solution is independent of terms first order in the perturbation parameter, and relax the assumptions needed for the local existence theorem of perturbation solutions.Perturbation, matrix calculus, DSGE, solution methods, Bézout theorem; Sylvester equations

Greedy low-rank algorithm for spatial connectome regression

Recovering brain connectivity from tract tracing data is an important computational problem in the neurosciences. Mesoscopic connectome reconstruction was previously formulated as a structured matrix regression problem (Harris et al., 2016), but existing techniques do not scale to the whole-brain setting. The corresponding matrix equation is challenging to solve due to large scale, ill-conditioning, and a general form that lacks a convergent splitting. We propose a greedy low-rank algorithm for connectome reconstruction problem in very high dimensions. The algorithm approximates the solution by a sequence of rank-one updates which exploit the sparse and positive definite problem structure. This algorithm was described previously (Kressner and Sirkovi\'c, 2015) but never implemented for this connectome problem, leading to a number of challenges. We have had to design judicious stopping criteria and employ efficient solvers for the three main sub-problems of the algorithm, including an efficient GPU implementation that alleviates the main bottleneck for large datasets. The performance of the method is evaluated on three examples: an artificial "toy" dataset and two whole-cortex instances using data from the Allen Mouse Brain Connectivity Atlas. We find that the method is significantly faster than previous methods and that moderate ranks offer good approximation. This speedup allows for the estimation of increasingly large-scale connectomes across taxa as these data become available from tracing experiments. The data and code are available online

Solving DSGE Models with a Nonlinear Moving Average

We introduce a nonlinear infinite moving average as an alternative to the standard state-space policy function for solving nonlinear DSGE models. Perturbation of the nonlinear moving average policy function provides a direct mapping from a history of innovations to endogenous variables, decomposes the contributions from individual orders of uncertainty and nonlinearity, and enables familiar impulse response analysis in nonlinear settings. When the linear approximation is saddle stable and free of unit roots, higher order terms are likewise saddle stable and first order corrections for uncertainty are zero. We derive the third order approximation explicitly and examine the accuracy of the method using Euler equation tests.Perturbation, nonlinear impulse response, DSGE, solution methods

Preconditioned low-rank Riemannian optimization for linear systems with tensor product structure

The numerical solution of partial differential equations on high-dimensional domains gives rise to computationally challenging linear systems. When using standard discretization techniques, the size of the linear system grows exponentially with the number of dimensions, making the use of classic iterative solvers infeasible. During the last few years, low-rank tensor approaches have been developed that allow to mitigate this curse of dimensionality by exploiting the underlying structure of the linear operator. In this work, we focus on tensors represented in the Tucker and tensor train formats. We propose two preconditioned gradient methods on the corresponding low-rank tensor manifolds: A Riemannian version of the preconditioned Richardson method as well as an approximate Newton scheme based on the Riemannian Hessian. For the latter, considerable attention is given to the efficient solution of the resulting Newton equation. In numerical experiments, we compare the efficiency of our Riemannian algorithms with other established tensor-based approaches such as a truncated preconditioned Richardson method and the alternating linear scheme. The results show that our approximate Riemannian Newton scheme is significantly faster in cases when the application of the linear operator is expensive.Comment: 24 pages, 8 figure

Preconditioning techniques for generalized Sylvester matrix equations

Sylvester matrix equations are ubiquitous in scientific computing. However, few solution techniques exist for their generalized multiterm version, as they now arise in an increasingly large number of applications. In this work, we consider algebraic parameter-free preconditioning techniques for the iterative solution of generalized multiterm Sylvester equations. They consist in constructing low Kronecker rank approximations of either the operator itself or its inverse. While the former requires solving standard Sylvester equations in each iteration, the latter only requires matrix-matrix multiplications, which are highly optimized on modern computer architectures. Moreover, low Kronecker rank approximate inverses can be easily combined with sparse approximate inverse techniques, thereby enhancing their performance with little or no damage to their effectiveness.Comment: 26 pages, 3 figures, 2 tables. Submitted manuscrip

Efficient Approaches for Enclosing the United Solution Set of the Interval Generalized Sylvester Matrix Equation

In this work, we investigate the interval generalized Sylvester matrix equation ${\bf{A}}X{\bf{B}}+{\bf{C}}X{\bf{D}}={\bf{F}}$ and develop some techniques for obtaining outer estimations for the so-called united solution set of this interval system. First, we propose a modified variant of the Krawczyk operator which causes reducing computational complexity to cubic, compared to Kronecker product form. We then propose an iterative technique for enclosing the solution set. These approaches are based on spectral decompositions of the midpoints of ${\bf{A}}$, ${\bf{B}}$, ${\bf{C}}$ and ${\bf{D}}$ and in both of them we suppose that the midpoints of ${\bf{A}}$ and ${\bf{C}}$ are simultaneously diagonalizable as well as for the midpoints of the matrices ${\bf{B}}$ and ${\bf{D}}$. Some numerical experiments are given to illustrate the performance of the proposed methods

Krylov subspace techniques for model reduction and the solution of linear matrix equations

This thesis focuses on the model reduction of linear systems and the solution of large scale linear matrix equations using computationally efficient Krylov subspace techniques. Most approaches for model reduction involve the computation and factorization of large matrices. However Krylov subspace techniques have the advantage that they involve only matrix-vector multiplications in the large dimension, which makes them a better choice for model reduction of large scale systems. The standard Arnoldi/Lanczos algorithms are well-used Krylov techniques that compute orthogonal bases to Krylov subspaces and, by using a projection process on to the Krylov subspace, produce a reduced order model that interpolates the actual system and its derivatives at infinity. An extension is the rational Arnoldi/Lanczos algorithm which computes orthogonal bases to the union of Krylov subspaces and results in a reduced order model that interpolates the actual system and its derivatives at a predefined set of interpolation points. This thesis concentrates on the rational Krylov method for model reduction. In the rational Krylov method an important issue is the selection of interpolation points for which various techniques are available in the literature with different selection criteria. One of these techniques selects the interpolation points such that the approximation satisfies the necessary conditions for H2 optimal approximation. However it is possible to have more than one approximation for which the necessary optimality conditions are satisfied. In this thesis, some conditions on the interpolation points are derived, that enable us to compute all approximations that satisfy the necessary optimality conditions and hence identify the global minimizer to the H2 optimal model reduction problem. It is shown that for an H2 optimal approximation that interpolates at m interpolation points, the interpolation points are the simultaneous solution of m multivariate polynomial equations in m unknowns. This condition reduces to the computation of zeros of a linear system, for a first order approximation. In case of second order approximation the condition is to compute the simultaneous solution of two bivariate polynomial equations. These two cases are analyzed in detail and it is shown that a global minimizer to the H2 optimal model reduction problem can be identified. Furthermore, a computationally efficient iterative algorithm is also proposed for the H2 optimal model reduction problem that converges to a local minimizer. In addition to the effect of interpolation points on the accuracy of the rational interpolating approximation, an ordinary choice of interpolation points may result in a reduced order model that loses the useful properties such as stability, passivity, minimum-phase and bounded real character as well as structure of the actual system. Recently in the literature it is shown that the rational interpolating approximations can be parameterized in terms of a free low dimensional parameter in order to preserve the stability of the actual system in the reduced order approximation. This idea is extended in this thesis to preserve other properties and combinations of them. Also the concept of parameterization is applied to the minimal residual method, two-sided rational Arnoldi method and H2 optimal approximation in order to improve the accuracy of the interpolating approximation. The rational Krylov method has also been used in the literature to compute low rank approximate solutions of the Sylvester and Lyapunov equations, which are useful for model reduction. The approach involves the computation of two set of basis vectors in which each vector is orthogonalized with all previous vectors. This orthogonalization becomes computationally expensive and requires high storage capacity as the number of basis vectors increases. In this thesis, a restart scheme is proposed which restarts without requiring that the new vectors are orthogonal to the previous vectors. Instead, a set of two new orthogonal basis vectors are computed. This reduces the computational burden of orthogonalization and the requirement of storage capacity. It is shown that in case of Lyapunov equations, the approximate solution obtained through the restart scheme approaches monotonically to the actual solution