One drawback to FunctionInterpolation on multivariate functions is that (apparently) it does no recursive subdivision to make a more accurate interpolation. (You can manually increase the density of the uniform grid with InterpolationPoints.)
Another approach would be to use the FEM functionality and ElementMeshInterpolation to construct the InterpolatingFunction. With a suitable "MeshRefinementFunction", one can adaptively refine the mesh to improve the accuracy of the interpolation as needed.
Here is an example using the OP's test function.
First a mesh refinement function is constructed that tells ToElementMesh to refine the mesh where the error is too great (see below for further discussion of the error). We also use an area criterion to stop refinement when the area of the triangles get too small (which also keeps the interpolation from getting unnecessarily big). The tolerance levels for each can be adjusted to suit one's needs.
Needs["NDSolve`FEM`"]
ClearAll[meshRefine];
Block[{x, y, z, m},
With[{f = Function[{y, x}, Sin[Sqrt[y - x]]], (* function to interpolate *)
interp = Simplify@ Quiet@ InterpolatingPolynomial[ (* linear interpolation-- *)
{{{y[[1]], x[[1]]}, z[[1]]}, (* the symbols x, y, z, m *)
{{y[[2]], x[[2]]}, z[[2]]}, (* becomes variables in Compile *)
{{y[[3]], x[[3]]}, z[[3]]}},
{m[[1]], m[[2]]}]},
meshRefine = Compile[{{vertices, _Real, 2}, {area, _Real, 0}},
Block[{x, y, z, m, f1, f2},
y = vertices[[All, 1]];
x = vertices[[All, 2]];
If[y[[3]] >= x[[3]] && y[[2]] >= x[[2]] && y[[1]] >= x[[1]],
z = f[y, x];
m = Mean[vertices]; (* use centroid as test point for error estimate *)
f1 = interp; (* use linear interpolation to estimate error *)
f2 = f[m[[1]], m[[2]]]; (* compare with function at centroid *)
Abs[f1 - f2] > 1*^-3 && (* precision/accuracy measure for the interpolation *)
area > 2*^-4, (* limit on subdivision of triangles *)
False]
]]]];
reg = ImplicitRegion[y >= x, {{y, 0, 10}, {x, 0, 10}}];
emesh = ToElementMesh[reg,
"ImproveBoundaryPosition" -> False, (* else boundary points don't satisfy y >= x *)
"MeshRefinementFunction" -> meshRefine]
testfn = ElementMeshInterpolation[{emesh},
Function[{y, x}, Sin[Sqrt[y - x]]] @@@ emesh["Coordinates"]]
Let's examine the results. The refinement is greater where f is changing more rapidly (high second derivative):
Show[emesh["Wireframe"], Frame -> True]
Plot3D[testfn[y, x], {y, x} ∈ emesh]
The refinement is fairly consistent measured by the accuracy of the interpolation. It is slightly greater at the boundary y == x because the refinement is stopped by the area criterion. The mesh created by ToElementMesh is quadratic, but we used linear interpolation to set a criterion for the error. (I suppose that for a sufficiently fine mesh, if the linear interpolation error estimate is $\epsilon$, then the actual interpolation error will be of the order $C\,\epsilon^2$, for a suitable constant $C$ that probably depends on the function f. In any case, we used $\epsilon=10^{-3}$ as a cut-off an got an interpolation with an error around $10^{-5}$.)
Plot3D[testfn[y, x] - Sin[Sqrt[y - x]] // RealExponent,
{y, x} ∈ emesh, PlotRange -> {-10, 0}]
f = FunctionInterpolation[Sin[Sqrt[y - ty]], {y, 0, 10}, {t, 0, 1}]? $\endgroup$ – BlacKow Apr 5 '16 at 17:58f[y,x/y]$\endgroup$ – BlacKow Apr 5 '16 at 18:10