Consider a series of random matrices $X\_n\in\mathbb{R}^{n\times m}$ consisting of i.i.d. entries, each with zero mean and variance $1/m$, and let $a\_n,b\_n\in\mathbb{R}^{n\times1}$ be two deterministic (or random and independent on $X$) vectors, say with bounded norm. 


I want to find the structure of some "nice"/"simple" "limit" function, $f\_n$, of the following term
$$
a\_n^T\left(X\_nX\_n^T+I\_n\right)^{-1}b\_n-f\_n\to0
$$ 
almost surely, as $n,m\to\infty$ with fixed ratio. 


**EDIT**: Due to Ofer's answer and comments I will consider some specific choice of $a\_n$:
$$
\frac{1}{n}w\_n^TX\_n^T\left(X\_nX\_n^T+I\_n\right)^{-1}b\_n-f\_n\to0
$$ 
Since $w\_n^TX\_n^T = \sum\_{i=1}^mw\_{i}x\_i^T$ where $x\_i$ is $i$th raw of $X\_n$, we can write that
\begin{align}
\frac{1}{n}w\_n^TX\_n^T\left(X\_nX\_n^T+I\_n\right)^{-1}b\_n = \frac{1}{n}\sum\_{i=1}^mw\_{i}x\_i^T\left(X\_nX\_n^T+I\_n\right)^{-1}b\_n
\end{align}
We know that $X\_nX\_n^T = \sum\_{i=1}^mx\_ix\_i^T$. Let $\left[X\_nX\_n^T\right]\_i = X\_nX\_n^T-x\_ix\_i^T$. Thus,
\begin{align}
\frac{1}{n}w\_n^TX\_n^T\left(X\_nX\_n^T+I\_n\right)^{-1}b\_n &= \frac{1}{n}\sum\_{i=1}^mw\_{i}x\_i^T\left(X\_nX\_n^T+I\_n\right)^{-1}b\_n\\
&=\frac{1}{n}\sum\_{i=1}^mw\_{i}\frac{x\_i^T\left(\left[X\_nX\_n^T\right]\_i+I\_n\right)^{-1}b\_n}{1+x\_i^T\left(\left[X\_nX\_n^T\right]\_i+I\_n\right)^{-1}x\_i}
\end{align}
Now, since $x\_i$ is independent on $\left(\left[X\_nX\_n^T\right]\_i+I\_n\right)^{-1}$, we know that (a.s.)
$$
x\_i^T\left(\left[X\_nX\_n^T\right]\_i+I\_n\right)^{-1}x\_i-\int (1+x)^{-1} \rho(dx)\to0
$$
where $\rho$ is the limit density of eigenvalues of $XX^T$. The same is true for the numerator. So, the speculation is that $f\_n$ behaves like 
\begin{align}
f\_n &= \frac{1}{n}\sum\_{i=1}^mw\_{i}\frac{x\_i^Tb\_n\int (1+x)^{-1} \rho(dx)}{1+\int (1+x)^{-1} \rho(dx)}\\
&=\frac{\int (1+x)^{-1} \rho(dx)}{1+\int (1+x)^{-1} \rho(dx)}\frac{1}{n}w\_n^TX\_n^Tb\_n
\end{align}


**Update**: Numerical calculations suggests that the above "limit" is not true, although, I can't really say that I completely understand where the rub is.



Determining the asymptotic behavior of some function of random matrix
For $a\_n=b\_n$ with $\|a\_n\|=\|b\_n\|$,
the result is standard (and can be found for example in papers of Bai and Silverstein, usually as a technical lemma in the appendix...): 
let $\rho$ be the limit density of eigenvalues of $XX^T$ (the Pastur-Marchenko law). Then
the limit you seek is asymptotically the normalized trace of $(XX^T+I)^{-1}$, i.e., $A:=\int (1+x)^{-1} \rho(dx)$.


If $a\_n\neq b\_n$ it requires a bit more work, but the answer should be $A\cdot \langle a\_n,b\_n\rangle$. The reasoning is similar: the expression you have is $\sum \lambda\_i 
\alpha\_i \beta\_i$ where $\alpha\_i$ are the coefficients of $a\_n$ in the eigenbase of 
$XX^T$. This will concentrate near its mean (to prove that, use results on the eigenvectors as in Silverstein's paper from the mid 90's on eigenvectors of covariance matrices; everything is much simpler in the Gaussian case), which will give the expression I wrote.


