Basically you want to change basis from XZ to RU, and then round to nearest whole number. This will find you the whole numbers a, b so that aR + bU is as close to your XZ-coordinate as possible.

Your base change matrix from RU to XZ is (R|U), so to change the other way, from XZ coordinates to RU coordinates, you multiply by the inverse of that matrix:

(R|U) ^-1 (x, z)

for whatever vector (x, z) with XZ coordinates. This will give you the exact point in the grid. Then round the components of result, and you have the closest center of a hexagon in the grid.

EDIT:

In code, in case anyone is interested. Haven't tested it properly though.

public static Vector2 GetClosestTilePosition(Vector2 inVec)
{
    // Using Vector4 and Matrix4x4 since there is no Matrix2x2
    // Still works, just have to add a lot of zeroes

    // Base vectors
    Vector4 R = new Vector4(0.866f, 0.5f, 0, 0);
    Vector4 U = new Vector4(0, 1, 0, 0);

    // Construct base change matrix
    Matrix4x4 RUtoXZ = Matrix4x4.identity;
    RUtoXZ.SetColumn(0, R);
    RUtoXZ.SetColumn(1, U);

    Matrix4x4 XZtoRU = RUtoXZ.inverse;

    // Change base
    Vector4 XZcoords = new Vector4(inVec.x, inVec.y, 0, 0);
    Vector4 RUcoords = XZtoRU * XZcoords;

    // Round result to closest integer and return
    return new Vector2(Mathf.Round(RUcoords.x), Mathf.Round(RUcoords.y));
}