We propose a non-commutative algorithm for multiplying 2x2 matrices using 7 coefficient products. This algorithm reaches simultaneously a better accuracy in practice compared to previously known such fast algorithms, and a time complexity bound with the best currently known leading term (obtained via alternate basis sparsification). To build this algorithm, we consider matrix and tensor norms bounds governing the stability and accuracy of numerical matrix multiplication. First, we reduce those bounds by minimizing a growth factor along the unique orbit of Strassen's 2x2-matrix multiplication tensor decomposition. Second, we develop heuristics for minimizing the number of operations required to realize a given bilinear formula, while further improving its accuracy. Third, we perform an alternate basis sparsification that improves on the time complexity constant and mostly preserves the overall accuracy.

Krylov methods rely on iterated matrix-vector products $A^k u_j$ for an $n\times n$ matrix $A$ and vectors $u_1,\ldots,u_m$. The space spanned by all iterates $A^k u_j$ admits a particular basis -- the \emph{maximal Krylov basis} -- which consists of iterates of the first vector $u_1, Au_1, A^2u_1,\ldots$, until reaching linear dependency, then iterating similarly the subsequent vectors until a basis is obtained. Finding minimal polynomials and Frobenius normal forms is closely related to computing maximal Krylov bases. The fastest way to produce these bases was, until this paper, Keller-Gehrig's 1985 algorithm whose complexity bound $O(n^\omega \log(n))$ comes from repeated squarings of $A$ and logarithmically many Gaussian eliminations. Here $\omega>2$ is a feasible exponent for matrix multiplication over the base field. We present an algorithm computing the maximal Krylov basis in $O(n^\omega\log\log(n))$ field operations when $m \in O(n)$, and even $O(n^\omega)$ as soon as $m\in O(n/\log(n)^c)$ for some fixed real $c>0$. As a consequence, we show that the Frobenius normal form together with a transformation matrix can be computed deterministically in $O(n^\omega (\log\log(n))^2)$, and therefore matrix exponentiation~$A^k$ can be performed in the latter complexity if $\log(k) \in O(n^{\omega-1-\varepsilon})$ for some fixed $\varepsilon>0$. A key idea for these improvements is to rely on fast algorithms for $m\times m$ polynomial matrices of average degree $n/m$, involving high-order lifting and minimal kernel bases.

Quasiseparable matrices are a class of rank-structured matrices widely used in numerical linear algebra and of growing interest in computer algebra, with applications in e.g. the linearization of polynomial matrices. Various representation formats exist for these matrices that have rarely been compared.

Proofs of Retrievability are protocols which allow a Client to store data remotely and to efficiently ensure, via audits, that the entirety of that data is still intact. Dynamic Proofs of Retrievability (DPoR) also support efficient retrieval and update of any small portion of the data. We propose a novel protocol for arbitrary outsourced data storage that achieves both low remote storage size and audit complexity. A key ingredient, that can be also of intrinsic interest, reduces to efficiently evaluating a secret polynomial at given public points, when the (encrypted) polynomial is stored on an untrusted Server. The Server performs the evaluations and also returns associated certificates. A Client can check that the evaluations are correct using the certificates and some pre-computed keys, more efficiently than re-evaluating the polynomial. Our protocols support two important features: the polynomial itself can be encrypted on the Server, and it can be dynamically updated by changing individual coefficients cheaply without redoing the entire setup. Our methods rely on linearly homomorphic encryption and pairings, and our implementation shows good performance for polynomial evaluations with millions of coefficients, and efficient DPoR with terabytes of data. For instance, for a 1TB database, compared to the state of art, we can reduce the Client storage by 5000x, communication size by 20x, and client-side audit time by 2x, at the cost of one order of magnitude increase in server-side audit time.

Multiplicity code decoders are based on Hermite polynomial interpolation with error correction. In order to have a unique Hermite interpolant one assumes that the field of scalars has characteristic 0 or ≥ 𝓁 +1, where 𝓁 is the maximum order of the derivatives in the list of values of the polynomial and its derivatives which are interpolated. For scalar fields of characteristic 𝓁+1, the minimum number of values for interpolating a polynomial of degree ≤ D is D+1+2E(𝓁+1) when ≤ E of the values are erroneous. Here we give an error-correcting Hermite interpolation algorithm that requires fewer values, that is, that can tolerate more errors, assuming that the characteristic of the scalar field is either 0 or ≥ D+1. Our algorithm requires (𝓁+1)D + 1 - (𝓁+1)𝓁/2 + 2E values.

New algorithms are presented for computing annihilating polynomials of Toeplitz, Hankel, and more generally Toeplitz+Hankel-like matrices over a field. Our approach follows works on Coppersmith's block Wiedemann method with structured projections, which have been recently successfully applied for computing the bivariate resultant. A first baby steps/giant steps approach --directly derived using known techniques on structured matrices-- gives a randomized Monte Carlo algorithm for the minimal polynomial of an (n x n) Toeplitz or Hankel-like matrix of displacement rank α using(Õnw-c(w) Õ c(w)) arithmetic operations, where (w) is the exponent of matrix multiplication and (c(2.373) = 0.523) for the best known value of (w). For generic Toeplitz+Hankel-like matrices a second algorithm computes the characteristic polynomial; in particular, when the displacement rank is considered constant, its cost is (Õn2-1/w). Previous algorithms required (O(n2) operations while the exponents presented here are respectively less than 1.86 and 1.58 with the best known estimate for (w).

We present a non-commutative algorithm for the multiplication of a 2 × 2-block-matrix by its transpose using 5 block products (3 recursive calls and 2 general products) over C or any field of prime characteristic. We use geometric considerations on the space of bilinear forms describing 2 × 2 matrix products to obtain this algorithm and we show how to reduce the number of involved additions. The resulting algorithm for arbitrary dimensions is a reduction of multiplication of a matrix by its transpose to general matrix product, improving by a constant factor previously known reductions. Finally we propose schedules with low memory footprint that support a fast and memory efficient practical implementation over a prime field. To conclude, we show how to use our result in L · D · LT factorization.

Proof of Retrievability (PoR) and Provable Data Possession (PDP) schemes have been proposed to ensure the integrity of stored data on untrusted servers. A successful PoR audit ensures, with high probability, that every piece of stored data is recoverable by the server. Most PoR schemes proposed have focused on bandwidth and computation cost, but in some deployment scenarios the size of remote storage can be the most expensive factor. We propose a simple PoR scheme which is a variant on some existing PDP work. Compared to existing audit routines, the server computation cost and bandwidth are higher, but the server storage cost is minimal. Our preliminary work indicates that deploying this scheme may be less costly in commercial cloud settings, depending on the cost structure and frequency of audits.

We present new algorithms to detect and correct errors in the lower-upper factorization of a matrix, or the triangular linear system solution, over an arbitrary field. Our main algorithms do not require any additional information or encoding other than the original inputs and the erroneous output. Their running time is softly linear in the dimension times the number of errors when there are few errors, smoothly growing to the cost of fast matrix multiplication as the number of errors increases. We also present applications to general linear system solving.

We present a novel recursive algorithm for reducing a symmetric matrix to a triangular factorization which reveals the rank profile matrix. That is, the algorithm computes a factorization P TA P = L D L T where P is a permutation matrix, L is lower triangular with a unit diagonal and D is symmetric block diagonal with 1 x 1 and 2 x 2 antidiagonal blocks. This algorithm requires O(n2rømega-2) arithmetic operations, with n the dimension of the matrix, r its rank and ømega an admissible exponent for matrix multiplication. Furthermore, experimental results demonstrate that our algorithm has very good performance: its computational speed matches that of its numerical counterpart and is twice as fast as the unsymmetric exact Gaussian factorization. By adapting the pivoting strategy developed in the unsymmetric case, we show how to recover the rank profile matrix from the permutation matrix and the support of the block-diagonal matrix. We also note that there is an obstruction in characteristic 2 for revealing the rank profile matrix, which requires to relax the shape of the block diagonal by allowing the 2-dimensional blocks to have a non-zero bottom-right coefficient. This relaxed decomposition can then be transformed into a standard PLDLTP T decomposition at a negligible cost.

In this paper, we give novel certificates for triangular equivalence and rank profiles. These certificates enable to verify the row or column rank profiles or the whole rank profile matrix faster than recomputing them, with a negligible overall overhead. We first provide quadratic time and space non-interactive certificates saving the logarithmic factors of previously known ones. Then we propose interactive certificates for the same problems whose Monte Carlo verification complexity requires a small constant number of matrix-vector multiplications, a linear space, and a linear number of extra field operations. As an application we also give an interactive protocol, certifying the determinant of dense matrices, faster than the best previously known one.

Consider solving a black box linear system, A(u) x = b(u), where the entries are polynomials in u over a field K, and A(u) is full rank. The solution, x = 1/g(u) f(u), where g is always the least common monic denominator, can be found by evaluating the system at distinct points ξl in K. The solution can be recovered even if some evaluations are erroneous. In [Boyer and Kaltofen, Proc. SNC 2014] the problem is solved with an algorithm that generalizes Welch/Berlekamp decoding of an algebraic Reed-Solomon code. Their algorithm requires the sum of a degree bound for the numerators plus a degree bound for the denominator of the solution. It is possible that the degree bounds input to their algorithm grossly overestimate the actual degrees. We describe an algorithm that given the same inputs uses possibly fewer evaluations to compute the solution. We introduce a second count for the number of evaluations required to recover the solution based on work by Stanley Cabay. The Cabay count includes bounds for the highest degree polynomial in the coefficient matrix and right side vector, but does not require solution degree bounds. Instead our algorithm iterates until the Cabay termination criterion is reached. At this point our algorithm returns the solution. Assuming we have the actual degrees for all necessary input parameters, we give the criterion that determines when the Cabay count is fewer than the generalized Welch/Berlekamp count.

The class of quasiseparable matrices is defined by a pair of bounds, called the quasiseparable orders, on the ranks of the sub-matrices entirely located in their strictly lower and upper triangular parts. These arise naturally in applications, as e.g. the inverse of band matrices, and are widely used for they admit structured representations allowing to compute with them in time linear in the dimension. We show, in this paper, the connection between the notion of quasiseparability and the rank profile matrix invariant, presented in [Dumas & al. ISSAC'15]. This allows us to propose an algorithm computing the quasiseparable orders (rL,rU) in time O{n2sω-2} where s=max(rL,rU) and ω the exponent of matrix multiplication. We then present two new structured representations, a binary tree of PLUQ decompositions, and the Bruhat, using respectively O{ns log n/s and O{ns} field elements instead of O{ns2} for the classical generator and O{ns log n} for the hierarchically semiseparable representations. We present algorithms computing these representations in time O{n2sω-2}. These representations allow a matrix-vector product in time linear in the size of their representation. Lastly we show how to multiply two such structured matrices in time O{n2sω-2}.

The row (resp. column) rank profile of a matrix describes the stair case shape of its row (resp. column) echelon form. In an ISSAC'13 paper, we proposed a recursive Gaussian elimination that can compute simultaneously the row and column rank profiles of a matrix, as well as those of all of its leading sub-matrices, in the same time as state of the art Gaussian elimination algorithms. Here we first study the conditions making a Gaussian elimination algorithm reveal this information. We propose the definition of a new matrix invariant, the rank profile matrix, summarizing all information on the row and column rank profiles of all the leading sub-matrices. We also explore the conditions for a Gaussian elimination algorithm to compute all or part of this invariant, through the corresponding PLUQ decomposition. As a consequence, we show that the classical iterative CUP decomposition algorithm can actually be adapted to compute the rank profile matrix. Used, in a Crout variant, as a base-case to our ISSAC'13 implementation, it delivers a significant improvement in efficiency. Second, the row (resp. column) echelon form of a matrix are usually computed via different dedicated triangular decompositions. We show here that, from some PLUQ decompositions, it is possible to recover the row and column echelon forms of a matrix and of any of its leading sub-matrices thanks to an elementary post-processing algorithm.

Exact linear algebra is a core component of many symbolic and algebraic computations, as it often delivers competitive theoretical complexities and also better harnesses the efficiency of modern computing infrastructures. In this tutorial we will present an overview on the recent advances in exact linear algebra algorithmic and implementation techniques, and highlight the few key ideas that have proven successful in their design. As an illustration, we will study in more details the computation of some matrix normal forms over a finite field or the ring of polynomials, specific to computer algebra.

We present algorithms performing sparse univariate polynomial interpolation with errors in the evaluations of the polynomial. Based on the initial work by Comer, Kaltofen and Pernet [Proc. ISSAC 2012], we define the sparse polynomial interpolation codes and state that their minimal distance is precisely the code-word length divided by twice the sparsity. At ISSAC 2012, we have given a decoding algorithm for as much as half the minimal distance and a list decoding algorithm up to the minimal distance.

Gaussian elimination with full pivoting generates a PLUQ matrix decomposition. Depending on the strategy used in the search for pivots, the permutation matrices can reveal some information about the row or the column rank profiles of the matrix. We propose a new pivoting strategy that makes it possible to recover at the same time both row and column rank profiles of the input matrix and of any of its leading sub-matrices. We propose a rank-sensitive and quad-recursive algorithm that computes the latter PLUQ triangular decomposition of an m x n matrix of rank r in O(mnrω-2) field operations, with ω the exponent of matrix multiplication. Compared to the LEU decomposition by Malashonock, sharing a similar recursive structure, its time complexity is rank sensitive and has a lower leading constant. Over a word size finite field, this algorithm also improves the practical efficiency of previously known implementations.

We propose algorithms performing sparse interpolation with errors, based on Prony's--Ben-Or's & Tiwari's algorithm, using a Berlekamp/Massey algorithm with early termination. First, we present an algorithm that can recover a t-sparse polynomial f from a sequence of values, where some of the values are wrong, spoiled by either random or misleading errors. Our algorithm requires bounds T ≥ t and E ≥ e, where e is the number of evaluation errors. It interpolates f(ωi) for i = 1,..., 2T(E + 1), where ω is a field element at which each non-zero term evaluates distinctly.

We study algorithm based fault tolerance techniques for supporting malicious errors in distributed computations based on Chinese remainder theorem. The description holds for both computations with integers or with polynomials over a field. It unifies the approaches of redundant residue number systems and redundant polynomial systems through the Reed Solomon decoding algorithm proposed by Gao. We propose several variations on the application of the extended Euclid algorithm, where the error correction rate is adaptive. Several improvements are studied, including the use of various criterions for the termination of the Euclidean Algorithm, and an acceleration using the Half-GCD techniques. When there is some redundancy in the input, a gap in the quotient sequence is stated at the step matching the error correction, which enables early termination parallel computations. Experiments are shown to compare these approaches.

International audience

We present algorithms and heuristics to compute the characteristic polynomial of a matrix given its minimal polynomial. The matrix is represented as a black-box, i.e., by a function to compute its matrix-vector product. The methods apply to matrices either over the integers or over a large enough finite field. Experiments show that these methods perform efficiently in practice. Combined in an adaptive strategy, these algorithms reach significant speedups in practice for some integer matrices arising in an application from graph theory.

In the past two decades, some major efforts have been made to reduce exact (e.g. integer, rational, polynomial) linear algebra problems to matrix multiplication in order to provide algorithms with optimal asymptotic complexity. To provide efficient implementations of such algorithms one need to be careful with the underlying arithmetic. It is well known that modular techniques such as the Chinese remainder algorithm or the p -adic lifting allow very good practical performance, especially when word size arithmetic is used. Therefore, finite field arithmetic becomes an important core for efficient exact linear algebra libraries. In this article, we study high performance implementations of basic linear algebra routines over word size prime fields: especially matrix multiplication; our goal being to provide an exact alternate to the numerical BLAS library. We show that this is made possible by a careful combination of numerical computations and asymptotically faster algorithms. Our kernel has several symbolic linear algebra applications enabled by diverse matrix multiplication reductions: symbolic triangularization, system solving, determinant, and matrix inverse implementations are thus studied.

A new randomized algorithm is presented for computing the characteristic polynomial of an n x n matrix over a field. Over a suffciently large field the asymptotic expected complexity of the algorithm is O(nθ)field operations, improving by a factor of log n on the worst case complexity of Keller-Gehrig's algorithm [11].

We deal with the computation of the characteristic polynomial of dense matrices over word size finite fields and over the integers. We first present two algorithms for finite fields: one is based on Krylov iterates and Gaussian elimination. We compare it to an improvement of the second algorithm of Keller-Gehrig. Then we show that a generalization of Keller-Gehrig's third algorithm could improve both complexity and computational time. We use these results as a basis for the computation of the characteristic polynomial of integer matrices. We first use early termination and Chinese remaindering for dense matrices. Then a probabilistic approach, based on integer minimal polynomial and Hensel factorization, is particularly well suited to sparse and/or structured matrices.

The FFLAS project has established that exact matrix multiplication over finite fields can be performed at the speed of the highly optimized numerical BLAS routines. Since many algorithms have been reduced to use matrix multiplication in order to be able to prove an optimal theoretical complexity, this paper shows that those optimal complexity algorithms, such as LSP factorization, rank determinant and inverse computation can also be the most efficient.

In this paper we study different implementations of finite field arithmetic, essential foundation of computer algebra. We focus on Galois fields of word size cardinality at most, with any characteristic. Classical representations as machine integers, floating point numbers, polynomials and Zech logarithms are compared. Furthermore, very efficient implementations of finite field dot products, matrix-vector products and matrix-matrix products (namely the symbolic equivalent of level 1, 2 and 3 BLAS) are presented. Our implementations have many symbolic linear algebra applications: symbolic triangularization, system solving, exact determinant computation, matrix normal form are such examples.