A recursive way of constructing preconditioning matrices for the stiffness matrix in the discretization of selfadjoint second order elliptic boundary value problems is proposed. It is based on a sequence of nested finite element spaces with the usual nodal basis functions. Using a node-ordering corresponding to the nested meshes, the finite element stiffness matrix is recursively split up into two-level block structures and is factored approximately in such a way that any successive Schur complement is replaced (approximated) by a matrix defined recursively and therefore only implicitly given. To solve a system with this matrix we need to perform a fixed number (\(\nu)\) of iterations on the preceding level using as an iteration matrix the preconditioning matrix already defined on that level. It is shown that by a proper choice of iteration parameters it suffices to use \(\nu >(1-\gamma^ 2)^{-1/2}\) iterations for the so constructed \(\nu\)-fold \(V\)-cycle (where \(\nu =2\) corresponds to a \(W\)-cycle) preconditioning matrices to be spectrally equivalent to the stiffness matrix. The conditions involve only the constant \(\gamma\) in the strengthened Cauchy-Bunyakovskii-Schwarz inequality for the corresponding two-level hierarchical basis function spaces and are therefore independent of the regularity of the solution for instance. If we use successive uniform refinements of the meshes the method is of optimal order of computational complexity, if \(\gamma^ 2<8/9\). Under reasonable assumptions on the finite element mesh, the condition numbers turn out to be so small that there are in practice few reasons to use an accelerated iterative method like the conjugate gradient method, for instance.