Efficiently solving many sets of linear equations without inversion or factorization Suppose I have the normal set of linear equations $Ax = b$. If I can store and manipulate $A$ I have a variety of techniques available to me such as inversion, factorization, or an iterative method. The first two are particularly useful if I want to solve not for just $x$ and $b$, but for a set of $x$ and $b$ vectors ($x_1, x_2,... x_n$ and $b_1, b_2,..., b_n$).
However, if $A$ is so large that I cannot explicitly store it and can only form $Ax$ products, inversion is certainly out and I am unaware of factorization techniques that only require $Ax$ products.
So my question is very simple: If I use an iterative method to solve for $x_1$ does this give me any information about the structure of $A$ that can be used to solve for $x_2,...$ more efficiently? If no, are there factorization techniques available that only require $Ax$ products?
A is both positive-definite and symmetric. I am currently solving for each $x_n$ via preconditioned conjugate gradient which is quite slow. Hopefully I have overlooking something obvious.
 A: We assume the classical requirements for success of a conjugate gradient method, i.e. $A$ is symmetric positive definite.  Use of a pre-conditioner with the modifications discussed below is not precluded.
A modified conjugate gradient method which deals with multiple right hand sides is found in the literature, called block conjugate gradient, going back to O'Leary (1980):

If there are $s$ systems to be solved, the block conjugate gradient algorithm will solve them in at most $\lceil n/s \rceil$ iterations and may involve less work than applying the conjugate gradient algorithm $s$ times.

The prospect for reduction in work in this original paper is tied to checking for rank deficiencies in the block residuals.  Since then various alternative approaches to checking for linear dependencies among the multiple residuals have been explored.
A survey of then current alternatives is presented by Dubrulle (2001) in Retooling the method of block conjugate gradients.
The most recent paper I spotted in this area is Block variants of the COCG and COCR methods for solving complex symmetric linear systems with multiple right-hand sides by Gu et al (2016):

In the present study, we establish two new block variants of the Conjugate
  Orthogonal Conjugate Gradient (COCG) and the Conjugate A-Orthogonal
  Conjugate Residual (COCR) Krylov subspace methods for solving complex symmetric
  linear systems with multiple right hand sides.

The ideas of block residuals can be exploited in more general Krylov space solvers.  Cf. Gutknecht's paper Block Krylov Space Methods for Linear Systems with Multiple Right-Hand Sides: An Introduction, 2006.

In a number of applications in scientific computing and engineering one has to solve huge sparse linear systems of equations with several right-hand sides that are given at once.

A: Conjugate gradient (and similar iterative methods) don't even know that there is a matrix involved, so it's hard to see how they could give us any information about its structure. 
Here's an idea (which I have never tried myself) ...
After you have solved a few different versions of the problem with different values of $b$, you have a few known values of the mapping $f(b) = A^{-1}b$. You could use these values to construct an approximation $g$ of $f$. Perhaps a polynomial or a spline approximation. Then, when you have to solve a problem with a new $b$, you use an iterative method, starting at $g(b)$.
Finding a good starting point is often the key to good speed and reliability in an iterative numerical algorithm. But, in this particular case, I'm not sure that a good starting point will help all that much. The problem is linear, so I would expect the iterations to converge very quickly, anyway. 
