A variational approach to unified microscopic description of normal and superfluid phases of a strongly interacting Bose system is proposed. We begin the formulation of an optimal theory within this approach through the diagrammatic analysis and synthesis of key distribution functions that characterize the spatial structure and the degree of coherence present in the two phases. The approach centers on functional minimization of the free energy corresponding to a suitable trial form for the many-body density matrix W(R, R') proportional to Phi(R) Q(R, R') Phi(R'), with the wave function Phi and incoherence factor Q chosen to incorporate the essential dynamical and statistical correlations. In earlier work addressing the normal phase, Phi was taken as a Jastrow product of two-body dynamical correlation factors and Phi was taken as a permanent of short-range two-body statistical bonds. A stratagem applied to the noninteracting Bose gas by Ziff, Uhlenbeck, and Kac is invoked to extend this ansatz to encompass both superfluid and normal phases, while introducing a variational parameter B that signals the presence of off-diagonal long-range order. The formal development proceeds from a generating Functional Lambda, defined by the logarithm of the normalization integral integral dR Phi(2)(R) Q(R, R). Construction of the Ursell-Mayer diagrammatic expansion of the generator Lambda is followed by renormalization of the statistical bond and of the parameter B. For B = 0, previous results for the normal phase are reproduced, whereas For B > 0, corresponding to the superfluid regime, a new class of anomalous contributions appears. Renormalized expansions for the pair distribution function g(r) and the cyclic distribution function G(cc)(r) are extracted from Lambda by functional differentiation. Standard diagrammatic techniques are adapted to obtain the appropriate hypernetted-chain equations for the evaluation of these spatial distribution functions. Corresponding results are presented for the internal energy. The quantity G(cc)(r) is found to develop long-range order in the condensed phase and therefore, assumes an incisive diagnostic role in the elucidation of the Bose-Einstein transition of the interacting system. A tentative connection of the microscopic description with the phenomenological two-fluid model is established in terms of a sum rule on ''normal'' and ''anomalous'' density components. Further work within this correlated density matrix approach will address the one-body density matrix, the entropy, and the Euler-Lagrange equations that lead to an optimal theory. (C) 1995 Academic Press, Inc.