Rapidly rotating, slightly nonaxisymmetric neutron stars emit nearly periodic gravitational waves (GWs), quite possibly at levels detectable by ground-based GW interferometers. We refer to these sources as ``GW pulsars.'' For any given sky position and frequency evolution, the $\mathcal{F}$-statistic is the maximum likelihood statistic for the detection of GW pulsars. However, in ``all-sky'' searches for previously unknown GW pulsars, it would be computationally intractable to calculate the (fully coherent) $\mathcal{F}$-statistic at every point of (a suitably fine) grid covering the parameter space: the number of grid points is many orders of magnitude too large for that. Therefore, in practice some nonoptimal detection statistic is used for all-sky searches. Here we introduce a ``phase-relaxed'' $\mathcal{F}$-statistic, which we denote ${\mathcal{F}}_{\mathrm{pr}}$, for incoherently combining the results of fully coherent searches over short time intervals. We estimate (very roughly) that for realistic searches, our ${\mathcal{F}}_{\mathrm{pr}}$ is $\ensuremath{\sim}10--15%$ more sensitive than the ``semicoherent'' $\mathcal{F}$-statistic that is currently used. Moreover, as a by-product of computing ${\mathcal{F}}_{\mathrm{pr}}$, one obtains a rough determination of the time-evolving phase offset between one's template and the true signal imbedded in the detector noise. Almost all the ingredients that go into calculating ${\mathcal{F}}_{\mathrm{pr}}$ are already implemented in the LIGO Algorithm Library, so we expect that relatively little additional effort would be required to develop a search code that uses ${\mathcal{F}}_{\mathrm{pr}}$.