Abstract: A numerically robust algorithm for computing ellipsoidal bounds on the state of a linear, discrete-time dynamic system was proposed. The algorithm employed ellipsoidal outer approximation of the feasible set assuming instantaneous process and observation noise vectors and the initial state to be bounded by known ellipsoids. The time and observation updates produced, respectively, the vector sum and intersection of two ellipsoids. Cholesky decomposition was used in the propagation of the shape-defining matrix of the ellipsoid to keep it positive definite in the presence of roundoff errors. Besides, a subminimal-volume ellipsoid was selected from a family of ellipsoids as the observation-updated ellipsoid to circumvent the complex optimization affected by ill-conditioned matrix inverse. Monte Carlo simulations on a digital computer were performed to compare the performance of the proposed algorithm with that of the optimal algorithm. Simulation results show that the proposed algorithm not only matches the performance of the optimal algorithm closely in terms of ellipsoid volumes and mean-square errors, but also is less vulnerable to roundoff errors. The proposed algorithm also features the capability to be realized on a parallel computer.