I agree with BrianBorchers' general approach. To summarize my discussion with him in the comment thread below his answer, I suggest some modifications to his approach:
To recap the gist of Brian's answer, the quickest way to calculate the trace of $A^{+}B$ for many values of $B$ is to precompute $A^{+}$, and then perform only the matrix-vector products needed for the diagonal elements of $A^{+}B$, and sum these diagonal elements. In the limit of large numbers of $B$, the expense of those matrix-vector products will dominate the cost of forming $A^{+}$. $A$ is also sparse, so you could exploit sparsity for speed, too. Memory becomes a concern when $A$ is large. You mentioned that $A$ is no larger than $5000 \times 10000$, so with double precision arithmetic, it could require up to 400 MB to store $A$. That's a lot, but storing it as a dense matrix is possible; also storing each $B$ as a dense matrix would not be possible.
Brian suggests calculating $A^{+}$ as $(A^{T}A)^{-1}A = A^{+}$, since $A$ has full column rank. We can rearrange this problem into $(A^{T}A)A^{+} = A$. This immediately suggests three approaches, all based on the normal equations of linear regression (i.e., least squares):
1) Use sparse Cholesky factorization of $A^{T}A$ to find a sparse matrix $R$ such that $R^{T}R = A^{T}A$; this step can be carried out using a call to CHOLMOD. Call $RA^{+} = X$. Solve the equation $R^{T}X = A$ for $X$ using a sparse LU factorization such as in UMFPACK; $R^{T}$ is lower triangular, so this solve should be quick, and only use symbolic factorization and some pivoting. Then solve $RA^{+} = X$ for $A^{+}$, again using a sparse LU factorization. $R$ is upper triangular, and it might be possible to reuse some of the information from the previous LU factorization in this solve.
2) Use a sparse QR factorization to find $R$ -- it's the same matrix as the one above. SuiteSparseQR is one implementation of sparse QR by Tim Davis, author of UMFPACK and CHOLMOD. This part saves you a matrix multiply (which doesn't really matter in the grand scheme of things), and you can either use the two sparse LU factorization approach above, or treat it as a least squares problem directly.
3) Use an SVD factorization of $A$, and the formula in this Wikipedia article.
Of these three approaches, the sparse QR decomposition and SVD approach are what I used in my thesis. In your case, I wouldn't necessarily recommend the SVD approach because I'm not aware of a good recent implementation of sparse SVD; if one is available, it's worth trying. Cholesky factorization is supposed to do a better job of exploiting sparsity than QR, and is probably fastest. QR is generally used in practice. SVD has the best numerical stability, followed by Cholesky and QR.
I assume in Brian's answer that when he suggests calculating $(A^{T}A)^{-1} = (R^{T}R)^{-1}$, he probably really means, "solve the system $(R^{T}R)A^{+} = A$". Solving a linear system by calculating the inverse explicitly and then multiplying is conventionally considered bad practice. Evidence suggests it's really not that bad for accuracy, but the inverse of a sparse matrix is frequently dense, so inverting the matrix is likely to require more storage and more operations than a sparse factorization. Furthermore, $R$ is better conditioned than $R^{T}R$, so recasting the linear system $(R^{T}R)A^{+} = A$ as two separate linear systems is at least as accurate as solving the single linear system; solving it as two separate systems is the conventional advice given in least squares problems.
Brian also rightly pointed out that there are approaches to calculate $A^{+}B$ directly using QR factorizations, and that these approaches would be unfavorable because they would require more computation and probably more memory.
