We study the problem of recovering an unknown signal $${\varvec{x}}$$
x
given measurements obtained from a generalized linear model with a Gaussian sensing matrix. Two popular solutions are based on a linear estimator $$\hat{\varvec{x}}^\mathrm{L}$$
x
^
L
and a spectral estimator $$\hat{\varvec{x}}^\mathrm{s}$$
x
^
s
. The former is a data-dependent linear combination of the columns of the measurement matrix, and its analysis is quite simple. The latter is the principal eigenvector of a data-dependent matrix, and a recent line of work has studied its performance. In this paper, we show how to optimally combine $$\hat{\varvec{x}}^\mathrm{L}$$
x
^
L
and $$\hat{\varvec{x}}^\mathrm{s}$$
x
^
s
. At the heart of our analysis is the exact characterization of the empirical joint distribution of $$({\varvec{x}}, \hat{\varvec{x}}^\mathrm{L}, \hat{\varvec{x}}^\mathrm{s})$$
(
x
,
x
^
L
,
x
^
s
)
in the high-dimensional limit. This allows us to compute the Bayes-optimal combination of $$\hat{\varvec{x}}^\mathrm{L}$$
x
^
L
and $$\hat{\varvec{x}}^\mathrm{s}$$
x
^
s
, given the limiting distribution of the signal $${\varvec{x}}$$
x
. When the distribution of the signal is Gaussian, then the Bayes-optimal combination has the form $$\theta \hat{\varvec{x}}^\mathrm{L}+\hat{\varvec{x}}^\mathrm{s}$$
θ
x
^
L
+
x
^
s
and we derive the optimal combination coefficient. In order to establish the limiting distribution of $$({\varvec{x}}, \hat{\varvec{x}}^\mathrm{L}, \hat{\varvec{x}}^\mathrm{s})$$
(
x
,
x
^
L
,
x
^
s
)
, we design and analyze an approximate message passing algorithm whose iterates give $$\hat{\varvec{x}}^\mathrm{L}$$
x
^
L
and approach $$\hat{\varvec{x}}^\mathrm{s}$$
x
^
s
. Numerical simulations demonstrate the improvement of the proposed combination with respect to the two methods considered separately.