  This work is on a fast and accurate reduced basis method for solving
discretized fractional elliptic partial differential equations (PDEs) of the
form $\mathcal{A}^su=f$ by rational approximation. A direct computation of the
action of such an approximation would require solving multiple (20$\sim$30)
large-scale sparse linear systems. Our method constructs the reduced basis
using the first few directions obtained from the preconditioned conjugate
gradient method applied to one of the linear systems. As shown in the theory
and experiments, only a small number of directions (5$\sim$10) are needed to
approximately solve all large-scale systems on the reduced basis subspace. This
reduces the computational cost dramatically because: (1) We only use one of the
large-scale problems to construct the basis; and (2) all large-scale problems
restricted to the subspace have much smaller sizes. We test our algorithms for
fractional PDEs on a 3d Euclidean domain, a 2d surface, and random
combinatorial graphs. We also use a novel approach to construct the rational
approximation for the fractional power function by the orthogonal greedy
algorithm (OGA).
