  In the least-squares fitting framework, the Vandermonde with Arnoldi (V+A)
method presented in [Brubeck, Nakatsukasa, and Trefethen, SIAM Review, 63
(2021), pp. 405-415] is an effective approach to compute a polynomial that
approximates an underlying univariate function f. Extensions of V+A include its
multivariate version and the univariate confluent V+A; the latter enables us to
use the information of the derivative of f in obtaining the approximation
polynomial. In this paper, we shall extend V+A further to the multivariate
confluent V+A. Besides the technical generalization of the univariate confluent
V+A, we also introduce a general and application-dependent G-orthogonalization
in the Arnoldi process. We shall demonstrate with several applications that, by
specifying an application-related G-inner product, the desired approximate
multivariate polynomial as well as its certain partial derivatives can be
computed accurately from a well-conditioned least-squares problem whose
coefficient matrix is orthonormal. The desired multivariate polynomial is
represented in a discrete G-orthogonal polynomials basis which admits an
explicit recurrence, and therefore, facilitates evaluating function values and
certain partial derivatives at new nodes efficiently. We demonstrate its
flexibility by applying it to solve the multivariate Hermite least-squares
problem and PDEs with various boundary conditions in irregular domains.
