In this manuscript, we propose a stable algorithm for computing the zeros of Althammer polynomials. These polynomials are orthogonal with respect to a Sobolev inner product, and are even if their degree is even, odd otherwise. Furthermore, their zeros are real, distinct, and located inside the interval (−1, 1). The Althammer polynomial p_n(x) of degree n satisfies a long recurrence relation, whose coefficients can be arranged into a Hessenberg matrix of order n, with eigenvalues equal to the zeros of the considered polynomial. Unfortunately, the eigenvalues of this Hessenberg matrix are very ill–conditioned, and standard balancing procedures do not improve their condition numbers. Here, we introduce a novel algorithm for computing the zeros of p_n(x), which first transforms the Hessenberg matrix into a similar symmetric tridiagonal one, i.e., a matrix whose eigenvalues are perfectly conditioned, and then computes the zeros of p_n(x) as the eigenvalues of the latter tridiagonal matrix. Moreover, we propose a second algorithm, faster but less accurate than the former one, which computes the zeros of p_n(x) as the eigenvalues of a truncated Hessenberg matrix, obtained by properly neglecting some diagonals in the upper part of the original matrix. The computational complexity of the proposed algorithms are, respectively, O( n^3/ 6 ), and O(l^2 n), with l<

In this work, we focus on the computation of the zeros of a monic Laguerre–Sobolev orthogonal polynomial of degree n. Taking into account the associated four–term recurrence relation, this problem can be formulated as a generalized eigenvalue problem, involving a lower bidiagonal matrix and a 2–banded lower Hessenberg matrix of order n. Unfortunately, the considered generalized eigenvalue problem is very ill–conditioned, and classical balancing procedures do not improve it. Therefore, customary techniques for solving the generalized eigenvalue problem, like the QZ method, yield unreliable results. Here, we propose a novel balancing procedure that drastically reduces the ill–conditioning of the eigenvalues of the involved matrix pencil. Moreover, we propose a fast and reliable algorithm, with O(n2) computational complexity and O(n) memory, exploiting the structure of the considered matrix pencil.

A fast and weakly stable method for computing the zeros of a particular class of hypergeometric polynomials is presented. The studied hypergeometric polynomials satisfy a higher order differential equation and generalize Laguerre polynomials. The theoretical study of the asymptotic distribution of the spectrum of these polynomials is an active research topic. In this article we do not contribute to the theory, but provide a practical method to contribute to further and better understanding of the asymptotic behavior. The polynomials under consideration fit into the class of Sobolev orthogonal polynomials, satisfying a four–term recurrence relation. This allows computing the roots via a generalized eigenvalue problem. After condition enhancing similarity transformations, the problem is transformed into the computation of the eigenvalues of a comrade matrix, which is a symmetric tridiagonal modified by a rank–one matrix. The eigenvalues are then retrieved by relying on an existing structured rank based fast algorithm. Numerical examples are reported studying the accuracy, stability and conforming the efficiency for various parameter settings of the proposed approach.

It is known that executing a perfect shifted QR step via the implicit QR algorithm may not result in a deflation of the perfect shift. Typically, several steps are required before deflation actually takes place. This deficiency can be remedied by determining the similarity transformation via the associated eigenvector. Similar techniques have been deduced for the QZ algorithm and for the rational QZ algorithm. In this paper we present a similar approach for executing a perfect shifted QZ step on a general rank structured pencil instead of a specific rank structured one, e.g., a Hessenberg--Hessenberg pencil. For this, we rely on the rank structures present in the transformed matrices. A theoretical framework is presented for dealing with general rank structured \rev{pencils} and deflating subspaces. We present the corresponding algorithm allowing} to deflate simultaneously a block of eigenvalues rather than a single one. We define the level-rho poles and show that these poles are maintained executing the deflating algorithm. Numerical experiments illustrate the robustness of the presented approach showing the importance of using the improved scaled residual approach.

In this paper we consider the computation of the modified moments for the system of Laguerre polynomials on the real semiaxis with the Hermite weight. These moments can be used for the computation of integrals with the Hermite weight on the real semiaxis via product rules. We propose a new computational method based on the construction of the null-space of a rectangular matrix derived from the three-term recurrence relation of the system of orthonormal Laguerre polynomials.It is shown that the proposed algorithm computes the modified moments with high relative accuracy and linear complexity.Numerical examples illustrate the effectiveness of the proposed method.

In this paper we analyze the stability of the problem of performing a rational QZ$step with a shift that is an eigenvalue of a given regular pencil H-lambda K in unreduced Hessenberg-Hessenberg form. In exact arithmetic, the backward rational QZ step moves the eigenvalue to the top of the pencil, while the rest of the pencil is maintained in Hessenberg-Hessenberg form, which then yields a deflation of the given shift. But in finite-precision the rational QZ step gets ``blurred'' and precludes the deflation of the given shift at the top of the pencil. In this paper we show that when we first compute the corresponding eigenvector to sufficient accuracy, then the rational QZ step can be constructed using this eigenvector, so that the exact deflation is also obtained in finite-precision.

The aim of this paper is to describe a Matlab package for computing the simultaneous Gaussian quadrature rules associated with a variety of multiple orthogonal polynomials. Multiple orthogonal polynomials can be considered as a generalization of classical orthogonal polynomials, satisfying orthogonality constraints with respect to different measures, with Moreover, they satisfy -term recurrence relations. In this manuscript, without loss of generality, is considered equal to The so-called simultaneous Gaussian quadrature rules associated with multiple orthogonal polynomials can be computed by solving a banded lower Hessenberg eigenvalue problem. Unfortunately, computing the eigendecomposition of such a matrix turns out to be strongly ill-conditioned and the Matlab function balance.m does not improve the condition of the eigenvalue problem. Therefore, most procedures for computing simultaneous Gaussian quadrature rules are implemented with variable precision arithmetic. Here, we propose a Matlab package that allows to reliably compute the simultaneous Gaussian quadrature rules in floating point arithmetic. It makes use of a variant of a new balancing procedure, recently developed by the authors of the present manuscript, that drastically reduces the condition of the Hessenberg eigenvalue problem.

In this paper, we derive a new method to compute the nodes and weights of simultaneous n-point Gaussian quadrature rules. The method is based on the eigendecomposition of the banded lower Hessenberg matrix that contains the coefficients of the recurrence relations for the corresponding multiple orthogonal polynomials. The novelty of the approach is that it uses the property of total nonnegativity of this matrix associated with the particular considered multiple orthogonal polynomials, in order to compute its eigenvalues and eigenvectors in a numerically stable manner. The overall complexity of the computation of all the nodes and weights is O(n^2).

The computation of n-point Gaussian quadrature rules for symmetric weight functions is considered in this paper. It is shown that the nodes and the weights of the Gaussian quadrature rule can be retrieved from the singular value decomposition of a bidiagonal matrix of size n/2. The proposed numerical method allows to compute the nodes with high relative accuracy and a computational complexity of O(n). We also describe an algorithm for computing the weights of a generic Gaussian quadrature rule with high relative accuracy. Numerical examples show the effectiveness of the proposed approach.

The aim of this work is to propose a fast and reliable algorithm for computingintegrals of the type$$\int_{-\infty}^{\infty} f(x) e^{\scriptstyle -x^2 -\frac{\scriptstyle 1}{\scriptstyle x^2}} dx,$$where $f(x)$ is a sufficiently smooth function, in floating point arithmetic.The algorithm is based on a product integration rule, whose rate of convergencedepends only on the regularity of $f$, since the coefficients of the rule are ``exactly'' computed by means of suitable recurrence relations here derived.We prove stability and convergence in the space of locally continuous functions on $\RR$ equipped with weighted uniform norm.By extensive numerical tests, the accuracy of the proposed product rule is compared with that of the Gauss--Hermite quadrature formula w.r.t. the function $f(x) e^{-\frac{\scriptstyle 1}{\scriptstyle x^2}}$. The numerical results confirm the effectiveness of the method, supporting the proven theoretical estimates.

This paper presents a sparse factorization for the delay Vandermonde matrix (DVM) and a faster, exact, radix-2, and completely recursive DVM algorithm to realize millimeter wave beamformers in wireless communication networks. The proposed algorithm will reduce the complexity of $N$-beam wideband beamformers from $\mathcal{O}(N^2)$ to $\mathcal{O}(N {\rm\: log\:} N)$. The scaled DVM algorithm is at least 97$\%$ faster than the brute-force scale DVM by a vector product. The signal flow graphs of the scaled DVM algorithm are shown to elaborate the simplicity of the proposed algorithm. The proposed lower complexity DVM algorithm can be used to design simple signal flow graph and realize in very large scale integrated circuit architecture with the significant reduction of chip area and power consumption. Moreover, the realization of the faster DVM algorithm through analog integrated circuits will be addressed . Finally, the proposed DVM algorithm will be utilized to obtain a low-complexity approximate transform for beamforming

The computation of the eigenvalue decomposition of matrices is one of the most investigated problems in numerical linear algebra. In particular, real nonsymmetric tridiagonal eigenvalue problems arise in a variety of applications. In this paper the problem of computing an eigenvector corresponding to a known eigenvalue of a real nonsymmetric tridiagonal matrix is considered, developing an algorithm that combines part of a QR sweep and part of a QL sweep, both with the shift equal to the known eigenvalue. The numerical tests show the reliability of the proposed method.

In this paper we revisit the problem of performing a QZ step with a so-called perfect shift, which is an exact eigenvalue of a given regular pencil lambdaB-A in unreduced Hessenberg triangular form. In exact arithmetic, the QZ step moves that eigenvalue to the bottom of the pencil, while the rest of the pencil is maintained in Hessenberg triangular form, which then yields a deflation of the given eigenvalue. But in finite precision the QZ step gets blurred and precludes the deflation of the given eigenvalue. In this paper we show that when we first compute the corresponding eigenvector to sufficient accuracy, then the QZ step can be constructed using this eigenvector, so that the deflation is also obtained in finite precision. An important application of this technique is the computation of the index of a system of differential algebraic equations, since an exact deflation of the infinite eigenvalues is needed to impose correctly the algebraic constraints of such differential equations.

In this paper we describe how to swap two 2 × 2 blocks in a real Schur form and a generalized real Schur form. We pay special attention to the numerical stability of the method. We also illustrate the stability of our approach by a series of numerical tests.

The computation of the eigenvalue decomposition of symmetric matrices is one of the most investigated problems in numerical linear algebra. For a matrix of moderate size, the customary procedure is to reduce it to a symmetric tridiagonal one by means of an orthogonal similarity transformation and then compute the eigendecomposition of the tridiagonal matrix. Recently, Malyshev and Dhillon have proposed an algorithm for deflating the tridiagonal matrix, once an eigenvalue has been computed. Starting from the aforementioned algorithm, in this manuscript we develop a procedure for computing an eigenvector of a symmetric tridiagonal matrix, once its associate eigenvalue is known. We illustrate the behavior of the proposed method with a number of numerical examples.

The paper proposes a decision support system (DSS) for the supply chain of packaged fresh and highly perishable products. The DSS combines a unique tool for sales forecasting with order planning which includes an individual model selection system equipped with ARIMA, ARIMAX and transfer function forecasting model families, the latter two accounting for the impact of prices. Forecasting model parameters are chosen via two alternative tuning algorithms: a two-step statistical analysis, and a sequential parameter optimisation framework for automatic parameter tuning. The DSS selects the model to apply according to user-defined performance criteria. Then, it considers sales forecasting as a proxy of expected demand and uses it as input for a multi-objective optimisation algorithm that defines a set of non-dominated order proposals with respect to outdating, shortage, freshness of products and residual stock. A set of real data and a benchmark - based on the methods already in use - are employed to evaluate the performance of the proposed DSS. The analysis of different configurations shows that the DSS is suitable for the problem under investigation; in particular, the DSS ensures acceptable forecasting errors and proper computational effort, providing order plans with associated satisfactory performances.

The generalized Schur algorithm is a powerful tool allowing to compute classical decompositions of matrices, such as the QR and LU factorizations. When applied to matrices with particular structures, the generalized Schur algorithm computes these factorizations with a complexity of one order of magnitude less than that of classical algorithms based on Householder or elementary transformations. In this manuscript, we describe the main features of the generalized Schur algorithm. We show that it helps to prove some theoretical properties of the R factor of the QR factorization of some structured matrices, such as symmetric positive definite Toeplitz and Sylvester matrices, that can hardly be proven using classical linear algebra tools. Moreover, we propose a fast implementation of the generalized Schur algorithm for computing the rank of Sylvester matrices, arising in a number of applications. Finally, we propose a generalized Schur based algorithm for computing the null-space of polynomial matrices.

We address the problem of forecasting sales for fresh and highly perishable products, in the general context of supply chain management. The forecasting activity refers to the single item in a given store and started from a pre-processing phase for data analysis and normalization. Then data was used as input for a forecasting algorithm designed to be user interactive. We implemented three forecasting methods: ARIMA, ARIMAX and transfer function models. The exogenous components of the forecasting models took the impact of prices into account. The best configuration of these models is dynamically chosen via two alternative methods: (i) a two-step procedure, based on properly selected statistical indicators, (ii) a Sequential Parameter Optimization approach for automatic parameter tuning. The user or the decision maker at the store level should not be exposed to the complexity of the forecasting system which - for this reason - is designed to adaptively select the best model configuration at every forecast session, to be used for each item/store combination. A set of real data based on 19 small and medium sized stores and 156 fresh products was employed to evaluate both quality of forecasting results and their effects on the order planning activity, where sales forecasting is considered as a proxy of the expected demand. Some examples are reported and discussed. Our results confirm that there is no 'one-size-fits-all' forecasting model, whose performance strictly depends on the specific characteristics of the underlying data. This supports the adoption of a data-driven tool to automate the dynamic selection of the most appropriate forecasting model. (C) 2017 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights reserved.

In this paper we revisit the problem of performing a QR-step on an unreduced Hessenberg matrix H when we know an "exact" eigenvalue ?0 of H. Under exact arithmetic, this eigenvalue will appear on diagonal of the transformed Hessenberg matrix H~ and will be decoupled from the remaining part of the Hessenberg matrix, thus resulting in a deflation. But it is well known that in finite precision arithmetic the so-called perfect shift can get blurred and that the eigenvalue ?0 can then not be deflated and/or is perturbed significantly. In this paper, we develop a new strategy for computing such a QR step so that the deflation is almost always successful. We also show how to extend this technique to double QR-steps with complex conjugate shifts.

An algorithm for computing the antitriangular factorization of symmetric matrices, relying only on orthogonal transformations, was recently proposed. The computed antitriangular form straightforwardly reveals the inertia of the matrix. A block version of the latter algorithm was described in a different paper, where it was noticed that the algorithm sometimes fails to compute the correct inertia of the matrix.In this paper we analyze a possible cause of the failure of detecting the inertia and propose a procedure to recover it. Furthermore, we propose a different algorithm to compute the antitriangular factorization of a symmetric matrix that handles most of the singularities of the matrix at the very end of the algorithm.Numerical results are also given showing the reliability of the proposed algorithm.