Let @math be a linear matrix, or pencil, generated by given symmetric matrices @math of size @math with rational entries. The set of real vectors x such that the pencil is positive semidefinite is a convex semi-algebraic set called spectrahedron, described by a linear matrix inequality (LMI). We design an exact algorithm that, up to genericity assumptions on the input matrices, computes an exact algebraic representation of at least one point in the spectrahedron, or decides that it is empty. The algorithm does not assume the existence of an interior point, and the computed point minimizes the rank of the pencil on the spectrahedron. The degree @math of the algebraic representation of the point coincides experimentally with the algebraic degree of a generic semidefinite program associated to the pencil. We provide explicit bounds for the complexity of our algorithm, proving that the maximum number of arithmetic operations that are performed is essentially quadratic in a multilinear B 'ezout bound of @math . When @math (resp. @math ) is fixed, such a bound, and hence the complexity, is polynomial in @math (resp. @math ). We conclude by providing results of experiments showing practical improvements with respect to state-of-the-art computer algebra algorithms.
Hence, in the case of symmetric matrices, these results cannot be used anymore. Because of the structure of the matrix, the system defining the incidence variety involves too many equations; some of them being redundant. Hence, these redundancies need to be eliminated to characterize critical points on incidence varieties in a convenient way. In the case of Hankel matrices, the special structure of their kernel provides an efficient way to do that. This case study is done in @cite_36 . Yet, the problem of eliminating these redundancies remained unsolved in the general symmetric case and this is what we do in which is the starting point of the design of our dedicated algorithm.
@cite_36
2170042437 1988452081 2061781302 1522470213
Let H0, …, H n be m x m matrices with entries in Q and Hankel structure, i.e. constant skew diagonals. We consider the linear Hankel matrix H(x) = H0+x1H_1+…+xnHn and the problem of computing sample points in each connected component of the real algebraic set defined by the rank constraint rank (H(x))≤ r, for a given integer r ≤ m-1. Computing sample points in real algebraic sets defined by rank defects in linear matrices is a general problem that finds applications in many areas such as control theory, computational geometry, optimization, etc. Moreover, Hankel matrices appear in many areas of engineering sciences. Also, since Hankel matrices are symmetric, any algorithmic development for this problem can be seen as a first step towards a dedicated exact algorithm for solving semi-definite programming problems, i.e. linear matrix inequalities. Under some genericity assumptions on the input (such as smoothness of an incidence variety), we design a probabilistic algorithm for tackling this problem. It is an adaptation of the so-called critical point method that takes advantage of the special structure of the problem. Its complexity reflects this: it is essentially quadratic in specific degree bounds on an incidence variety. We report on practical experiments and analyze how the algorithm takes advantage of this special structure. A first implementation outperforms existing implementations for computing sample points in general real algebraic sets: it tackles examples that are out of reach of the state-of-the-art. The generalized eigenvalue problem @math with H a Hankel matrix and @math the corresponding shifted Hankel matrix occurs in number of applications such as the reconstruction of the shape of a polygon from its moments, the determination of abscissa of quadrature formulas, of poles of Pade approximants, or of the unknown powers of a sparse black box polynomial in computer algebra. In many of these applications, the entries of the Hankel matrix are only known up to a certain precision. We study the sensitivity of the nonlinear application mapping the vector of Hankel entries to its generalized eigenvalues. A basic tool in this study is a result on the condition number of Vandermonde matrices with not necessarily real abscissas which are possibly row-scaled. The goal of this paper is to present an efficient implementation of an explicit matrix inversion of general square matrices on multicore computer architecture. The inversion procedure is split into four steps: 1) computing the LU factorization, 2) inverting the upper triangular U factor, 3) solving a linear system, whose solution yields inverse of the original matrix and 4) applying backward column pivoting on the inverted matrix. Using a tile data layout, which represents the matrix in the system memory with an optimized cache-aware format, the computation of the four steps is decomposed into computational tasks. A directed acyclic graph is generated on the fly which represents the program data flow. Its nodes represent tasks and edges the data dependencies between them. Previous implementations of matrix inversions, available in the state-of-the-art numerical libraries, are suffer from unnecessary synchronization points, which are non-existent in our implementation in order to fully exploit the parallelism of the underlying hardware. Our algorithmic approach allows to remove these bottlenecks and to execute the tasks with loose synchronization. A runtime environment system called QUARK is necessary to dynamically schedule our numerical kernels on the available processing units. The reported results from our LU-based matrix inversion implementation significantly outperform the state-of-the-art numerical libraries such as LAPACK (5x), MKL (5x) and ScaLAPACK (2.5x) on a contemporary AMD platform with four sockets and the total of 48 cores for a matrix of size 24000. A power consumption analysis shows that our high performance implementation is also energy efficient and substantially consumes less power than its competitors. We describe new data structures for full and packed storage of dense symmetric triangular arrays that generalize both. Using the new data structures one is led to several new algorithms that save "half" the storage for symmetric matrices and outperform the current blocked based level 3 algorithms in LAPACK. We concentrate on the simplest forms of the new algorithms and show they are a direct generalization of LINPACK. This means that level 3 BLAS's are not required to obtain level 3 performance. The replacement for Level 3 BLAS are so-called kernel routines, see [1], and on IBM platforms they are producible from simple textbook type codes, by the XLF Fortran compiler. In the sequel I will label these "vanilla" codes. On Power3 with a peak performance of 800 MFlops, the results for Cholesky factorization at order n ? 200 is over 720 MFlops and then reaches 735 MFlops at n = 400. Using conventional full format LAPACK DPOTRF with ESSL BLAS's one first gets to 600 MFlops at n ? 600 and only reaches a peak of 620 MFlops. For this result we used simple square blocked full matrix data formats where the blocks themselves are stored in column major (Fortran) order or row major (C) order. The simple algorithms of LU factorization with partial pivoting for this new data format is a direct generalization of LINPACK algorithm DGEFA. Again, no conventional level 3 BLAS's are required; the replacements are again so-called kernel routines. Programming for squared blocked full matrix format can be accomplished in standard Fortran through the use of three and four dimensional arrays. Thus, no new compiler support is necessary. Also we mention that other more complicated algorithms are possible; e.g., recursive ones. The recursive algorithms are also easily programmed via the use of tables that address where the blocks are stored in the two dimensional recursive block array. Finally, we describe block hybrid formats. Doing so allows one to use no additional storage over conventional (full and packed) matrix storage. This means the new algorithms are completely portable.