  We propose a way to construct a thermal pure quantum matrix product state
(TPQ-MPS) that can simulate finite temperature quantum many-body systems with a
minimal numerical cost comparable to the matrix product algorithm for the
ground state. The MPS was originally designed for the wave function with
area-law entanglement. However, by attaching the auxiliary sites to the edges
of the random matrix product state, we find that the degree of entanglement is
automatically tuned so as to recover the volume law of the entanglement entropy
that characterizes the TPQ state. The finite temperature physical quantities of
the transverse Ising and the spin-1/2 Heisenberg chains evaluated by a TPQ-MPS
show excellent agreement even for bond dimension $\sim 10$-$20$ with those of
the exact results.
