scispace - formally typeset
Search or ask a question

Showing papers in "SIAM Journal on Matrix Analysis and Applications in 2012"


Journal ArticleDOI
TL;DR: This paper proposes that the random variation is best described via a Poisson distribution, which better describes the zeros observed in the data as compared to the typical assumption of a Gaussian distribution, and presents a new algorithm for Poisson tensor factorization called CANDECOMP--PARAFAC alternating Poisson regression (CP-APR), based on a majorization-minimization approach.
Abstract: Tensors have found application in a variety of fields, ranging from chemometrics to signal processing and beyond. In this paper, we consider the problem of multilinear modeling of sparse count data. Our goal is to develop a descriptive tensor factorization model of such data, along with appropriate algorithms and theory. To do so, we propose that the random variation is best described via a Poisson distribution, which better describes the zeros observed in the data as compared to the typical assumption of a Gaussian distribution. Under a Poisson assumption, we fit a model to observed data using the negative log-likelihood score. We present a new algorithm for Poisson tensor factorization called CANDECOMP--PARAFAC alternating Poisson regression (CP-APR) that is based on a majorization-minimization approach. It can be shown that CP-APR is a generalization of the Lee--Seung multiplicative updates. We show how to prevent the algorithm from converging to non-KKT points and prove convergence of CP-APR under mil...

250 citations


Journal ArticleDOI
TL;DR: A local convergence theorem for calculating canonical low-rank tensor approximations (PARAFAC, CANDECOMP) by the alternating least squares algorithm is established and the main assumption is that the Hessian matrix of the problem is positive definite modulo the scaling indeterminacy.
Abstract: A local convergence theorem for calculating canonical low-rank tensor approximations (PARAFAC, CANDECOMP) by the alternating least squares algorithm is established. The main assumption is that the Hessian matrix of the problem is positive definite modulo the scaling indeterminacy. A discussion, whether this is realistic, and numerical illustrations are included. Also regularization is addressed.

212 citations


Journal ArticleDOI
TL;DR: The problem of optimal model order reduction of bilinear control systems with respect to the generalization of the well-known ${\cal H}_2$-norm for linear systems is discussed.
Abstract: In this paper, we will discuss the problem of optimal model order reduction of bilinear control systems with respect to the generalization of the well-known ${\cal H}_2$-norm for linear systems. We...

162 citations


Journal ArticleDOI
TL;DR: An inductive method is introduced for the study of the uniqueness of decompositions of tensors, by means of Tensors of rank $1, based on the geometric notion of weak defectivity, which proves that the decomposition is unique for three-dimensional tensors of type a,b,c.
Abstract: We introduce an inductive method for the study of the uniqueness of decompositions of tensors, by means of tensors of rank $1$. The method is based on the geometric notion of weak defectivity. For three-dimensional tensors of type $(a,b,c)$, $a\leq b\leq c$, our method proves that the decomposition is unique (i.e., $k$-identifiability holds) for general tensors of rank $k$, as soon as $k\leq (a+1)(b+1)/16$. This improves considerably the known range for identifiability. The method applies also to tensor of higher dimension. For tensors of small size, we give a complete list of situations where identifiability does not hold. Among them, there are $4\times 4\times 4$ tensors of rank $6$, an interesting case because of its connection with the study of DNA strings.

161 citations


Journal ArticleDOI
TL;DR: This article motivates, derive, and test effective preconditioners to be used with the Minres algorithm for solving a number of saddle point systems which arise in PDE-constrained optimization problems.
Abstract: In this article, we motivate, derive, and test effective preconditioners to be used with the Minres algorithm for solving a number of saddle point systems which arise in PDE-constrained optimization problems. We consider the distributed control problem involving the heat equation and the Neumann boundary control problem involving Poisson's equation and the heat equation. Crucial to the effectiveness of our preconditioners in each case is an effective approximation of the Schur complement of the matrix system. In each case, we state the problem being solved, propose the preconditioning approach, prove relevant eigenvalue bounds, and provide numerical results which demonstrate that our solvers are effective for a wide range of regularization parameter values, as well as mesh sizes and time-steps.

115 citations


Journal ArticleDOI
TL;DR: A positive lower bound for the best rank-$1$ approximation ratio of a symmetric tensor is given and a higher order polynomial spherical optimization problem can be reformulated as a multilinear spherical optimized problem.
Abstract: In this paper, we show that for a symmetric tensor, its best symmetric rank-$1$ approximation is its best rank-$1$ approximation. Based on this result, a positive lower bound for the best rank-$1$ approximation ratio of a symmetric tensor is given. Furthermore, a higher order polynomial spherical optimization problem can be reformulated as a multilinear spherical optimization problem. Then, we present a modified power algorithm for solving the homogeneous polynomial spherical optimization problem. Numerical results are presented, illustrating the effectiveness of the proposed algorithm.

111 citations


Journal ArticleDOI
TL;DR: A superfast solver for Toeplitz linear systems based on rank structured matrix methods and randomized sampling is proposed, which uses displacement equations to transform a ToePlitz matrix.
Abstract: We propose a superfast solver for Toeplitz linear systems based on rank structured matrix methods and randomized sampling. The solver uses displacement equations to transform a Toeplitz matrix $T$ ...

109 citations


Journal ArticleDOI
TL;DR: This paper focuses on the construction of explicit low-rank representations of matrices in the tensor train (TT) and quantized tensorTrain (QTT) formats which have been proposed recently for the low-parametric structured representation of large-scale tensors.
Abstract: We focus on the construction of explicit low-rank representations of matrices in the tensor train (TT) and quantized tensor train (QTT) formats which have been proposed recently for the low-parametric structured representation of large-scale tensors. The matrices under consideration are discretizations of the Laplace operator in a hypercube (in one and many dimensions) and their inverses (in one dimension). For these matrices we derive explicit and exact QTT representations of low QTT ranks independent of the numbers $d$ and $n^{2}$ of dimensions and entries in each dimension. This implies that for the matrices considered the storage cost is $\mathcal{O}\left(d\,\log n\right)$, i.e., logarithmic with respect to the total number $n^{2d}$ of entries. The same applies to the computational complexity of the QTT-structured operations with these matrices, which now depends on the QTT ranks of the other operands. The general result of the paper is the notation and technique we introduce in order to examine the Q...

107 citations


Journal ArticleDOI
TL;DR: Orthogonality-constrained versions of the CPD methods based on simultaneous matrix diagonalization and alternating least squares are presented and a simple proof of the existence of the optimal low-rank approximation of a tensor in the case that a factor matrix is columnwise orthonormal is given.
Abstract: Canonical polyadic decomposition (CPD) of a higher-order tensor is an important tool in mathematical engineering. In many applications at least one of the matrix factors is constrained to be columnwise orthonormal. We first derive a relaxed condition that guarantees uniqueness of the CPD under this constraint. Second, we give a simple proof of the existence of the optimal low-rank approximation of a tensor in the case that a factor matrix is columnwise orthonormal. Third, we derive numerical algorithms for the computation of the constrained CPD. In particular, orthogonality-constrained versions of the CPD methods based on simultaneous matrix diagonalization and alternating least squares are presented. Numerical experiments are reported.

73 citations


Journal ArticleDOI
TL;DR: It is shown that the minimal indices of the pencil can be determined with only the Wong sequences and that the Kronecker canonical form is a simple corollary of this result; hence, in passing, this work provides a new proof for the Kr onecker canonicalform.
Abstract: We study singular matrix pencils and show that the so-called Wong sequences yield a quasi-Kronecker form. This form decouples the matrix pencil into an underdetermined part, a regular part, and an overdetermined part. This decoupling is sufficient to fully characterize the solution behavior of the differential-algebraic equations associated with the matrix pencil. Furthermore, we show that the minimal indices of the pencil can be determined with only the Wong sequences and that the Kronecker canonical form is a simple corollary of our result; hence, in passing, we also provide a new proof for the Kronecker canonical form. The results are illustrated with an example given by a simple electrical circuit.

72 citations


Journal ArticleDOI
TL;DR: This paper establishes concise notation that is suitable for the analysis and development of block Tensor algorithms, proves several useful block tensor identities, and makes precise the notion of a block tensorsor unfolding.
Abstract: Within the field of numerical multilinear algebra, block tensors are increasingly important. Accordingly, it is appropriate to develop an infrastructure that supports reasoning about block tensor computation. In this paper we establish concise notation that is suitable for the analysis and development of block tensor algorithms, prove several useful block tensor identities, and make precise the notion of a block tensor unfolding.

Journal ArticleDOI
TL;DR: Hierarchical structured matrices have been widely used in fast solutions of integral equations, PDEs, structured matrix (such as Toeplitz) problems, companion eigenproblems, etc.
Abstract: In recent years, hierarchical structured matrices have been widely used in fast solutions of integral equations, PDEs, structured matrix (such as Toeplitz) problems, companion eigenproblems, etc. I...

Journal ArticleDOI
TL;DR: This work assumes that the real matrix of the LCP is an $H-matrix and solves it by using a new method, the scaled extrapolated block modulus algorithm, as well as an improved version of the very recently introduced modulus-based matrix splitting modified AOR iteration method.
Abstract: The numerous applications of the linear complementarity problem (LCP) in, e.g., the solution of linear and convex quadratic programming, free boundary value problems of fluid mechanics, and moving boundary value problems of economics make its efficient numerical solution a very imperative and interesting area of research. For the solution of the LCP, many iterative methods have been proposed, especially, when the matrix of the problem is a real positive definite or an $H_{+}$-matrix. In this work we assume that the real matrix of the LCP is an $H_{+}$-matrix and solve it by using a new method, the scaled extrapolated block modulus algorithm, as well as an improved version of the very recently introduced modulus-based matrix splitting modified AOR iteration method. As is shown by numerical examples, the two new methods are very effective and competitive with each other. (A corrected PDF is attached to this article.)

Journal ArticleDOI
TL;DR: It is demonstrated is that all three methods are capable of delivering minimal nonnegative solutions with entrywise relative accuracies as warranted by the defining coefficient matrices of a MARE.
Abstract: A new doubling algorithm—the alternating-directional doubling algorithm (ADDA)—is developed for computing the unique minimal nonnegative solution of an $M$-matrix algebraic Riccati equation (MARE). It is argued by both theoretical analysis and numerical experiments that ADDA is always faster than two existing doubling algorithms: SDA of Guo, Lin, and Xu (Numer. Math., 103 (2006), pp. 393-412) and SDA-ss of Bini, Meini, and Poloni (Numer. Math., 116 (2010), pp. 553-578) for the same purpose. Also demonstrated is that all three methods are capable of delivering minimal nonnegative solutions with entrywise relative accuracies as warranted by the defining coefficient matrices of a MARE. The three doubling algorithms, differing only in their initial setups, correspond to three special cases of the general bilinear (also called Mobius) transformation. It is explained that ADDA is the best among all possible doubling algorithms resulted from all bilinear transformations.

Journal ArticleDOI
TL;DR: Two theoretical results are presented that mirror the well-known trace minimization principle and Cauchy's interlacing inequalities for the symmetric eigenvalue problem.
Abstract: We present two theoretical results for the linear response eigenvalue problem. The first result is a minimization principle for the sum of the smallest eigenvalues with the positive sign. The second result is Cauchy-like interlacing inequalities. Although the linear response eigenvalue problem is a nonsymmetric eigenvalue problem, these results mirror the well-known trace minimization principle and Cauchy's interlacing inequalities for the symmetric eigenvalue problem.

Journal ArticleDOI
TL;DR: It is proved that if there is uniqueness in one mode, then the initial CP model can be uniquely decomposed in a sum of lower-rank tensors for which identifiability can be independently assessed.
Abstract: In this paper, three sufficient conditions are derived for the three-way CANDECOMP/PARAFAC (CP) model, which ensure uniqueness in one of the three modes (“uni-mode-uniqueness”). Based on these conditions, a partial uniqueness condition is proposed which allows collinear loadings in only one mode. We prove that if there is uniqueness in one mode, then the initial CP model can be uniquely decomposed in a sum of lower-rank tensors for which identifiability can be independently assessed. This condition is simpler and easier to check than other similar conditions existing in the specialized literature. These theoretical results are illustrated by numerical examples.

Journal ArticleDOI
TL;DR: This work presents how to obtain explicit description of AE parametric solution sets by combining a modified Fourier--Motzkin type elimination of existentially quantified parameters with the elimination of the universally quantification parameters.
Abstract: Consider linear systems whose input data are linear functions of uncertain parameters varying within given intervals. We are interested in an explicit description of the so-called AE parametric solution sets (where all universally quantified parameters precede all existentially quantified ones) by a set of inequalities not involving the parameters. This work presents how to obtain explicit description of AE parametric solution sets by combining a modified Fourier--Motzkin type elimination of existentially quantified parameters with the elimination of the universally quantified parameters. Some necessary (and sufficient) conditions for existence of nonempty AE parametric solution sets are discussed, as well as some properties of the parametric AE solution sets, e.g., shape of the solution set and some inclusion relations. Explicit descriptions of particular classes of AE parametric solution sets (tolerable, controllable, any two-dimensional) are given. Numerical examples illustrate the solution sets and th...

Journal ArticleDOI
TL;DR: This paper studies complete positivity of matrices whose underlying graph possesses a specific sparsity pattern, for example, being acyclic or circular, where the underlying graph of a symmetric matrix of order $n$ is defined to be a graph with $ n$ vertices and an edge between two vertices if the corresponding entry in the matrix is nonzero.
Abstract: A matrix $X$ is called completely positive if it allows a factorization $X = \sum_{b\in \mathcal{B}} bb{^\mathsf{T}}$ with nonnegative vectors $b$. These matrices are of interest in optimization, as it has been found that several combinatorial and quadratic problems can be formulated over the cone of completely positive matrices. The difficulty is that checking complete positivity is $\mathcal{NP}$-hard. Finding a factorization of a general completely positive matrix is also hard. In this paper we study complete positivity of matrices whose underlying graph possesses a specific sparsity pattern, for example, being acyclic or circular, where the underlying graph of a symmetric matrix of order $n$ is defined to be a graph with $n$ vertices and an edge between two vertices if the corresponding entry in the matrix is nonzero. The types of matrices that we analyze include tridiagonal matrices as an example. We show that in these cases checking complete positivity can be done in linear-time. A factorization of ...

Journal ArticleDOI
TL;DR: In this article, the facial structure of the set of correlation matrices, a convex set also known as the elliptope, is determined. But the facial structures are not usually considered together.
Abstract: In this paper we establish links between, and new results for, three problems that are not usually considered together. The first is a matrix decomposition problem that arises in areas such as statistical modeling and signal processing: given a matrix $X$ formed as the sum of an unknown diagonal matrix and an unknown low-rank positive semidefinite matrix, decompose $X$ into these constituents. The second problem we consider is to determine the facial structure of the set of correlation matrices, a convex set also known as the elliptope. This convex body, and particularly its facial structure, plays a role in applications from combinatorial optimization to mathematical finance. The third problem is a basic geometric question: given points $v_1,v_2,\ldots,v_n\in \mathbb{R}^k$ (where $n > k$) determine whether there is a centered ellipsoid passing exactly through all the points. We show that in a precise sense these three problems are equivalent. Furthermore we establish a simple sufficient condition on a su...

Journal ArticleDOI
TL;DR: An approach to avoid diverging components for real $I\times J\times K$ arrays with $R\le\min(I,J,K)$.
Abstract: Fitting an R-component Candecomp/Parafac (CP) decomposition to a multiway array or higher-order tensor ${\cal Z}$ is equivalent to finding a best rank-R approximation of ${\cal Z}$. Such a best rank-R approximation may not exist due to the fact that the set of multiway arrays with rank at most R is not closed. In this case, trying to compute the approximation results in diverging CP components. We present an approach to avoid diverging components for real $I\times J\times K$ arrays with $R\le\min(I,J,K)$. We show that a CP decomposition $({\bf A},{\bf B},{\bf C})$ featuring diverging components can be rewritten as a decomposition in block terms, where each block term corresponds to a group of diverging components. Moreover, we show that if the diverging components occur in groups of two or three, then the limiting boundary point ${\cal X}$ (i.e., the limit of the sequence of CP updates) can be obtained by fitting an appropriate constrained Tucker3 model to ${\cal Z}$, using the block term decomposition of...

Journal ArticleDOI
TL;DR: Two new normwise bounds are obtained for discrete-time Markov chains, and it is shown how the bounds developed in this paper and one bound given in [E. Seneta, Adv. Appl. Probab., 20 (1988), pp. 228--230] can be extended to continuous-timeMarkov chains.
Abstract: In this paper, we are interested in investigating the perturbation bounds for the stationary distributions for discrete-time or continuous-time Markov chains on a countable state space. For discrete-time Markov chains, two new normwise bounds are obtained. The first bound is rather easy to obtain since the needed condition, equivalent to uniform ergodicity, is imposed on the transition matrix directly. The second bound, which holds for a general (possibly periodic) Markov chain, involves finding a drift function. This drift function is closely related to the mean first hitting times. Some $V$-normwise bounds are also derived based on the results in [N. V. Kartashov, J. Soviet Math., 34 (1986), pp. 1493--1498]. Moreover, we show how the bounds developed in this paper and one bound given in [E. Seneta, Adv. Appl. Probab., 20 (1988), pp. 228--230] can be extended to continuous-time Markov chains. Several examples are shown to illustrate our results or to compare our bounds with the known ones in the literature.

Journal ArticleDOI
TL;DR: It is shown that arbitrary convergence behavior of Ritz values is possible in the Arnoldi method, and two parametrizations of the class of matrices with initial Arnoldi vectors that generate prescribed Ritzvalues are given.
Abstract: We show that arbitrary convergence behavior of Ritz values is possible in the Arnoldi method, and we give two parametrizations of the class of matrices with initial Arnoldi vectors that generate prescribed Ritz values (in all iterations). The second parametrization enables us to prove that any GMRES residual norm history is possible with any prescribed Ritz values (in all iterations), provided that we treat the stagnation case appropriately.

Journal ArticleDOI
TL;DR: This paper develops computable a posteriori error estimates for the pointwise evaluation of linear functionals of a solution to a parameterized linear system of equations and defines an improved linear functional that converges at a much faster rate than the original linear functional given a pointwise convergence assumption.
Abstract: We develop computable a posteriori error estimates for the pointwise evaluation of linear functionals of a solution to a parameterized linear system of equations. These error estimates are based on a variational analysis applied to polynomial spectral methods for forward and adjoint problems. We also use this error estimate to define an improved linear functional and we prove that this improved functional converges at a much faster rate than the original linear functional given a pointwise convergence assumption on the forward and adjoint solutions. The advantage of this method is that we are able to use low order spectral representations for the forward and adjoint systems to cheaply produce linear functionals with the accuracy of a higher order spectral representation. The method presented in this paper also applies to the case where only the convergence of the spectral approximation to the adjoint solution is guaranteed. We present numerical examples showing that the error in this improved functional is often orders of magnitude smaller. We also demonstrate that in higher dimensions, the computational cost required to achieve a given accuracy is much lower using the improved linear functional.

Journal ArticleDOI
TL;DR: A preconditioning technique based on a differencing approach such that the preconditionsed covariance matrix has a bounded condition number independent of the size of the matrix for some important process classes is discussed.
Abstract: In many statistical applications one must solve linear systems involving large, dense, and possibly irregularly structured covariance matrices. These matrices are often ill-conditioned; for example, the condition number increases at least linearly with respect to the size of the matrix when observations of a random process are obtained from a fixed domain. This paper discusses a preconditioning technique based on a differencing approach such that the preconditioned covariance matrix has a bounded condition number independent of the size of the matrix for some important process classes. When used in large scale simulations of random processes, significant improvement is observed for solving these linear systems with an iterative method.

Journal ArticleDOI
TL;DR: This work shows that a general iteration of X_{k+1} = f(X_k) for computing the unitary polar factor is backward stable under two conditions and proves the backward stability of the scaled Newton iteration under the assumption of mixed backward--forward stable manner.
Abstract: Among the many iterations available for computing the polar decomposition the most practically useful are the scaled Newton iteration and the recently proposed dynamically weighted Halley iteration. Effective ways to scale these and other iterations are known, but their numerical stability is much less well understood. In this work we show that a general iteration $X_{k+1} = f(X_k)$ for computing the unitary polar factor is backward stable under two conditions. The first condition requires that the iteration is implemented in a mixed backward--forward stable manner and the second requires that the mapping $f$ does not significantly decrease the size of any singular value relative to the largest singular value. Using this result we show that the dynamically weighted Halley iteration is backward stable when it is implemented using Householder QR factorization with column pivoting and either row pivoting or row sorting. We also prove the backward stability of the scaled Newton iteration under the assumption ...

Journal ArticleDOI
TL;DR: In this article, the authors developed a fast and accurate algorithm for computing con-eigenvalues and coneigenvectors of positive-definite Cauchy matrices, yielding even the tiniest con eigenvalues with high relative accuracy.
Abstract: The need to compute small con-eigenvalues and the associated con-eigenvectors of positive-definite Cauchy matrices naturally arises when constructing rational approximations with a (near) optimally small $L^{\infty}$ error. Specifically, given a rational function with $n$ poles in the unit disk, a rational approximation with $m\ll n$ poles in the unit disk may be obtained from the $m$th con-eigenvector of an $n\times n$ Cauchy matrix, where the associated con-eigenvalue $\lambda_{m}>0$ gives the approximation error in the $L^{\infty}$ norm. Unfortunately, standard algorithms do not accurately compute small con-eigenvalues (and the associated con-eigenvectors) and, in particular, yield few or no correct digits for con-eigenvalues smaller than the machine roundoff. We develop a fast and accurate algorithm for computing con-eigenvalues and con-eigenvectors of positive-definite Cauchy matrices, yielding even the tiniest con-eigenvalues with high relative accuracy. The algorithm computes the $m$th con-eigenval...

Journal ArticleDOI
TL;DR: In this article, a technique for using the evolution of the system to tell us about its initial structure, and then use this technique to develop an algorithm that takes the varied solutions from multiple clustering algorithms to arrive at a single data clustering solution.
Abstract: In 1961 Herbert Simon and Albert Ando [Econometrika, 29 (1961), pp. 111--138] published the theory behind the long-term behavior of a dynamical system that can be described by a nearly uncoupled matrix. Over the past fifty years this theory has been used in a variety of contexts, including queueing theory, brain organization, and ecology. In all of these applications, the structure of the system is known and the point of interest is the various stages the system passes through on its way to some long-term equilibrium. This paper looks at this problem from the other direction. That is, we develop a technique for using the evolution of the system to tell us about its initial structure, and then use this technique to develop an algorithm that takes the varied solutions from multiple data clustering algorithms to arrive at a single data clustering solution.

Journal ArticleDOI
TL;DR: This paper develops a strategy involving orthogonal transformations and approximations which avoids the explicit computation of the Schur complement in each factorization step, and utilizes a robustness technique so that an approximate generalized Cholesky factorization is guaranteed to exist.
Abstract: In this paper, we propose a robust Cholesky factorization method for symmetric positive definite (SPD), hierarchically semiseparable (HSS) matrices. Classical Cholesky factorizations and some semiseparable methods need to sequentially compute Schur complements. In contrast, we develop a strategy involving orthogonal transformations and approximations which avoids the explicit computation of the Schur complement in each factorization step. The overall factorization requires fewer floating point operations and has better data locality when compared to the recent HSS method in [J. Xia and M. Gu, SIAM J. Matrix Anal. Appl., 31 (2010), pp. 2899--2920]. Our strategy utilizes a robustness technique so that an approximate generalized Cholesky factorization is guaranteed to exist. We test three different methods for compressing the off-diagonal blocks in each iteration, i.e., rank-revealing QR, SVD, and SVD with random sampling. In our comparisons, we find that, with high probability, using SVD with random samplin...

Journal ArticleDOI
TL;DR: This paper presents multishift variants of Francis's implicitly shifted QR algorithm and provides a convergence theory that shows that the new algorithm performs nested subspace iterations on rational Krylov subspaces and confirms the validity of the theory.
Abstract: Recently a generalization of Francis's implicitly shifted QR algorithm was proposed, notably widening the class of matrices admitting low-cost implicit QR steps. This unifying framework covered the methods and theory for Hessenberg and inverse Hessenberg matrices and furnished also new, single-shifted, QR-type methods for, e.g., CMV matrices. Convergence of this approach was only suggested by numerical experiments. No theoretical claims supporting the results were presented. In this paper we present multishift variants of these new algorithms. We also provide a convergence theory that shows that the new algorithm performs nested subspace iterations on rational Krylov subspaces. Numerical experiments confirm the validity of the theory.

Journal ArticleDOI
TL;DR: It is shown that a simple doubling algorithm using this representation can reach full machine accuracy on a wide range of problems, obtaining invariant subspaces of the same quality as those computed by the state-of-the-art algorithms based on orthogonal transformations.
Abstract: We derive a new representation of Lagrangian subspaces in the form ${\rm Im} \Pi^T\big[\begin{smallmatrix}I \\ X\end{smallmatrix}\big]$, where $\Pi$ is a symplectic matrix which is the product of a permutation matrix and a real orthogonal diagonal matrix, and $X$ satisfies $\left\vert X_{ij}\right\vert \leq \begin{cases}1 & \text{if $i=j$,}\\ \sqrt{2} & \text{if $i eq j$.} \end{cases}$ This representation allows us to limit element growth in the context of doubling algorithms for the computation of Lagrangian subspaces and the solution of Riccati equations. It is shown that a simple doubling algorithm using this representation can reach full machine accuracy on a wide range of problems, obtaining invariant subspaces of the same quality as those computed by the state-of-the-art algorithms based on orthogonal transformations. The same idea carries over to representations of arbitrary subspaces and can be used for other types of structured pencils.