Will PCA always fit a model at least as well as NMF? If I perform PCA/NMF on a dataset, and then use the reduced models to reconstruct the original dataset, it seems to me that PCA should typically outperform NMF, simply due to the fact that NMF has the additional restriction of all positive matrix elements. So, are there any guarantees that PCA will outperform (or match) NMF (for the same number of components), perhaps as measured by least-squares or Poisson loss between the original and reconstructed data, or by BIC/AIC?
In Python, an example workflow is:
from sklearn.decomposition import NMF, PCA

X = <some data matrix of shape (n_samples, n_features)>

model_PCA = PCA(n_components=3)
W_PCA = model_PCA.fit_transform(X)
X_recon_PCA = model_PCA.inverse_transform(W_PCA)

model_NMF = NMF(n_components=3)
W_NMF = model_NMF.fit_transform(X)
X_recon_NMF = model_NMF.inverse_transform(W_NMF)

What can we say about the "quality" of X_recon_PCA vs X_recon_NMF compared to the true data matrix X?

Yes, by definition the PCA reconstruction will outperform any other linear combination of $k$ components in terms of RSS as it is designed to get the best $k$-rank approximation to the matrix A in terms of least-squares. To quote Pearson's (1901) On lines and planes of closest fit to systems of points in space: "the best fitting plane is perpendicular to the greatest axis of the ellipsoid of residuals and the minimum mean square residuals varies inversely as the length of this axis." In contrast to that, NNMF components are not constrained to be perpendicular so do not return the optimal reconstruction in terms of RSS.  Please note that they are extensions to a specific member of the exponential family of distributions (e.g.  Bernoulli, Poisson, etc.), Collins et al. (2001) A Generalization of Principal Components Analysis to the Exponential Family is an easy to follow and influential reference if you want to read on this further.
