$$ \begin{aligned} P(y \mid X) &= \int \mathcal{N}(y \mid X \omega, \sigma^2I)\mathcal{N}(\omega \mid 0, A)d\omega \\ &= \mathcal{N}(y \mid 0, \sigma^2I + XAX^t) \end{aligned} $$

$$ \begin{aligned} P(y \mid X) &= \int \mathcal{N}(y \mid X \omega, \sigma^2I)\mathcal{N}(\omega \mid 0, A)d\omega \\ &= \int \tfrac{1}{(2\pi)^{N/2}\sigma^N}\exp\left(\tfrac{-1}{2\sigma^2}\|y - X\omega\|^2 \right)\tfrac{1}{(2\pi)^{K/2}\sqrt{|A|}}\exp\left(\tfrac{-1}{2}\omega^t A^{-1} \omega \right)d\omega \\ &= \tfrac{1}{(2\pi)^{N/2}\sigma^N}\tfrac{1}{(2\pi)^{K/2}\sqrt{|A|}} \int \exp\left(\tfrac{-1}{2\sigma^2}\|y - X\omega\|^2 +\tfrac{-1}{2}\omega^t A^{-1} \omega \right)d\omega \\ &= \tfrac{1}{(2\pi)^{N/2}\sigma^N}\tfrac{1}{(2\pi)^{K/2}\sqrt{|A|}} \int \exp\left(\tfrac{-1}{2\sigma^2}\left(y^ty +\omega^tX^tX\omega - 2\omega^tX^ty \right) +\tfrac{-1}{2}\omega^t A^{-1} \omega \right)d\omega \\ &= \tfrac{1}{(2\pi)^{N/2}\sigma^N}\tfrac{1}{(2\pi)^{K/2}\sqrt{|A|}} \int \exp\left(\tfrac{-1}{2\sigma^2}\left(y^ty +\omega^tX^tX\omega - 2\omega^tX^ty +\color{green}{\sigma^2}\omega^t A^{-1} \omega \right)\right)d\omega \\ &= \tfrac{1}{(2\pi)^{N/2}\sigma^N}\tfrac{1}{(2\pi)^{K/2}\sqrt{|A|}} \int \exp\left(\tfrac{-1}{2\sigma^2}\left(y^ty +\omega^t\left(X^tX+\sigma^2A^{-1}\right)\omega - 2\omega^tX^ty \right)\right)d\omega \\ &= \tfrac{1}{(2\pi)^{N/2}\sigma^N}\tfrac{1}{(2\pi)^{K/2}\sqrt{|A|}}\exp\left(\tfrac{-y^ty}{2\sigma^2}\right) \int \exp\left(\tfrac{-1}{2\sigma^2}\left(\omega^t\left(X^tX+\sigma^2A^{-1}\right)\omega - 2\omega^tX^ty \right)\right)d\omega \quad \text{(1)} \end{aligned} $$

Now if we manage to write all the $\omega$ dependencies of the integral in a quadratic form we can use:

$$ \int \exp\left( \tfrac{-1}{2 \sigma^2} (\omega-\omega^*)^t(X^tX+\sigma^2 A^{-1})(\omega-\omega^*)\right)d\omega = \tfrac{(2\pi)^{K/2}\sigma^N}{\sqrt{|X^tX+\sigma^2 A^{-1}|}} $$

Hence if we expand this term we get:

$$ \begin{aligned} \tfrac{(2\pi)^{K/2}\sigma^N}{\sqrt{|X^tX+\sigma^2 A^{-1}|}} &= \int \exp\left(\tfrac{-1}{2 \sigma^2}(\omega-\omega^*)^t(X^tX+\sigma^2 A^{-1})(\omega-\omega^*)\right)d\omega \\ &= \int \exp\left(\tfrac{-1}{2 \sigma^2}\left(\omega(X^tX+\sigma^2 A^{-1})\omega + \omega^{*t}(X^tX+\sigma^2 A^{-1})\omega^* -2\omega^{t}(X^tX+\sigma^2 A^{-1})\omega^{*} \right)\right)d\omega \\ &= \exp\left(\tfrac{-\omega^{*t}(X^tX+\sigma^2 A^{-1})\omega^*}{2 \sigma^2}\right)\int \exp\left(\tfrac{-1}{2 \sigma^2}\left(\omega(X^tX+\sigma^2 A^{-1})\omega -2\omega^{t}(X^tX+\sigma^2 A^{-1})\omega^{*} \right)\right)d\omega \quad \text{(2)} \end{aligned} $$

Therefore if we want that both expressions ressemble we must set $\omega^*$ to: $$ \begin{aligned} \require{cancel} \cancel{-2}(X^tX +\sigma^2 A^{-1})\omega^* &= \cancel{-2}X^ty \\ \omega^* &= (X^tX +\sigma^2 A^{-1})^{-1}X^ty \end{aligned} $$

Which is the mean of the posterior and the MAP estimator!

If we plug this $\omega^*$ in the equation of above (2) we get:

$$ \begin{aligned} \tfrac{(2\pi)^{K/2}\sigma^N}{\sqrt{|X^tX+\sigma^2 A^{-1}|}} &= \exp\left(\frac{-y^tX(X^tX + \sigma^2 A^{-1})^{-1}X^ty}{2 \sigma^2}\right) \int \exp\left(\tfrac{-1}{2\sigma^2}\left(\omega(X^tX+\sigma^2 A^{-1})\omega -2\omega^{t}X^ty \right)\right)d\omega \\ \tfrac{(2\pi)^{K/2}\sigma^N}{\sqrt{|X^tX+\sigma^2 A^{-1}|}} \exp\left(\frac{y^tX(X^tX + \sigma^2 A^{-1})^{-1}X^ty}{2 \sigma^2}\right) &= \int \exp\left(\tfrac{-1}{2\sigma^2}\left(\omega(X^tX+\sigma^2 A^{-1})\omega -2\omega^{t}X^ty \right)\right)d\omega \quad \text{(3)} \end{aligned} $$

So now we can plug the right hand side of (3) into (1) to get: $$ \begin{aligned} P(y \mid X) &= \tfrac{1}{(2\pi)^{N/2}\cancel{\sigma^N}}\tfrac{1}{\cancel{(2\pi)^{K/2}}\sqrt{|A|}}\exp\left(\tfrac{-y^ty}{2\sigma^2}\right) \tfrac{\cancel{(2\pi)^{K/2}} \cancel{\sigma^N}}{\sqrt{|X^tX+\sigma^2 A^{-1}|}} \exp\left(\frac{y^tX(X^tX + \sigma^2 A^{-1})^{-1}X^ty}{2 \sigma^2}\right) \\ &= \tfrac{1}{(2\pi)^{N/2}\sqrt{|A||X^tX+\sigma^2 A^{-1}|}} \exp\left(\tfrac{-1}{2\sigma^2}\left(y^ty -y^tX(X^tX + \sigma^2 A^{-1})^{-1}X^ty\right)\right) \\ &= \tfrac{1}{(2\pi)^{N/2}\sqrt{|A||X^tX+\sigma^2 A^{-1}|}} \exp\left(\tfrac{-1}{2\sigma^2}y^t\left(I -X(X^tX + \sigma^2 A^{-1})^{-1}X^t\right)y\right) \\ &=^{\color{red}{\text{??}}} \mathcal{N}\left(y \mid 0, \sigma^2\left(I -X(X^tX + \sigma^2 A^{-1})^{-1}X^t\right)^{-1}\right) \end{aligned} $$

The matrix inversion lemma applied to the expression in the exponential leads to: $$ \begin{aligned} (\sigma^2 I + XAX^t)^{-1} &= \tfrac{1}{\sigma^2}\left(I-X\left(A^{-1} +X^tX\sigma^2\right)^{-1}X^t \sigma^{-2}\right)\\ &= \tfrac{1}{\sigma^2}\left(I-X\left(A^{-1}\sigma^2 +X^tX\right)^{-1}X^t \right)\\ \end{aligned} $$

Which completes the proof: $$ P(y \mid X ) = \mathcal{N}(y \mid 0, \sigma^2I + XAX^t) $$