In many engineering applications it is required to compute the dominant subspace of a matrix A of dimension m
× n, with m ≫ n. Often the matrix A is produced incrementally, so all the columns are not available simultaneously. This problem arises, e.g., in image processing, where each column of the matrix A represents an image of agiven sequence leading to a singular value decomposition based compression . Furthermore, the so called proper orthogonal decomposition approximation uses the left dominant subspace of a matrix A where a column consists of a time instance of the solution of an evolution equation, e.g., the ﬂow ﬁeld from a ﬂuid dynamics simulation. Since these ﬂow ﬁelds tend to be very large, only a small number can be stored eﬃciently during the simulation, and therefore an incremental approach is useful .
In this paper an algorithm for computing an approximation of the left dominant subspace of size k of A
∈ Rm×n, with k≪ m, n, is proposed requiring at each iteration O(mk + k2 ) ﬂoating point operations. Moreover, the proposed algorithm exhibits a lot of parallelism that can be exploited for a suitable implementation on a parallel computer.