Polynomial interpolation is well understood on the real line. In
multi-dimensional spaces, one often adopts a well established one-dimensional method
and fills up the space using certain tensor product rule. Examples
like this include full tensor construction and sparse grids construction.
This approach typically results in
fast growth of the total number of interpolation nodes and certain fixed geometrical structure of the
nodal sets. This imposes difficulties for practical applications,
where obtaining function values at a large number of nodes is
infeasible. Also, one often has function data from nodal locations that are not
by "mathematical design" and are ``unstructured''.

In this talk, we present a mathematical framework for conducting polynomial interpolation
in multiple dimensions using arbitrary set of unstructured nodes. The resulting method,
least orthogonal interpolation, is rigorous and has a straightforward numerical implementation.
It can faithfully interpolate any function data on any nodal sets, even on those that are considered
degenerate by the traditional methods. We also present a strategy to choose
``optimal'' nodes that result in robust
interpolation. The strategy is based on optimization of Lebesgue
function and has certain highly desirable mathematical properties.