  Vandermonde matrices are usually exponentially ill-conditioned and often
result in unstable approximations. In this paper, we introduce and analyze the
\textit{multivariate Vandermonde with Arnoldi (V+A) method}, which is based on
least-squares approximation together with a Stieltjes orthogonalization
process, for approximating continuous, multivariate functions on
$d$-dimensional irregular domains. The V+A method addresses the
ill-conditioning of the Vandermonde approximation by creating a set of discrete
orthogonal basis with respect to a discrete measure. The V+A method is simple
and general. It relies only on the sample points from the domain and requires
no prior knowledge of the domain. In this paper, we first analyze the sample
complexity of the V+A approximation. In particular, we show that, for a large
class of domains, the V+A method gives a well-conditioned and near-optimal
$N$-dimensional least-squares approximation using $M=\mathcal{O}(N^2)$
equispaced sample points or $M=\mathcal{O}(N^2\log N)$ random sample points,
independently of $d$. We also give a comprehensive analysis of the error
estimates and rate of convergence of the V+A approximation. Based on the
multivariate V+A approximation, we propose a new variant of the weighted V+A
least-squares algorithm that uses only $M=\mathcal{O}(N\log N)$ sample points
to give a near-optimal approximation. Our numerical results confirm that the
(weighted) V+A method gives a more accurate approximation than the standard
orthogonalization method for high-degree approximation using the Vandermonde
matrix.
