  A new efficient numerical algorithm for interacting fermion systems is
proposed and examined in detail. The ground state is expressed approximately by
a linear combination of numerically chosen basis states in a truncated Hilbert
space. Two procedures lead to a better approximation. The first is a numerical
renormalization, which optimizes the chosen basis and projects onto the ground
state within the fixed dimension, $L$, of the Hilbert space. The second is an
increase of the dimension of the truncated Hilbert space, which enables the
linear combination to converge to a better approximation. The extrapolation
$L\to\infty$ after the convergence removes the approximation error
systematically. This algorithm does not suffer from the negative sign problem
and can be applied to systems in any spatial dimension and arbitrary lattice
structure. The efficiency is tested and the implementation explained for
two-dimensional Hubbard models where Slater determinants are employed as chosen
basis. Our results with less than 400 chosen basis indicate good accuracy
within the errorbar of the best available results as those of the quantum Monte
Carlo for energy and other physical quantities.
