CODEWITHSUNDEEP
pythonscipyinterpolation
Adrian
Adrian

Reputation: 3288

Python Rbf gives singular matrix error with no duplicate coordinates, why?

Very similar to RBF interpolation fails: LinAlgError: singular matrix but I think the problem is different, as I have no duplicated coordinates.

Toy example:

import numpy as np
import scipy.interpolate as interp
coords = (np.array([-1, 0, 1]), np.array([-2, 0, 2]), np.array([-1, 0, 1]))
coords_mesh = np.meshgrid(*coords, indexing="ij")
fn_value = np.power(coords_mesh[0], 2) + coords_mesh[1]*coords_mesh[2]  # F(x, y, z)
coords_array = np.vstack([x.flatten() for x in coords_mesh]).T  # Columns are x, y, z
unique_coords_array = np.vstack({tuple(row) for row in coords_array})
unique_coords_array.shape == coords_array.shape  # True, i.e. no duplicate coords
my_grid_interp = interp.RegularGridInterpolator(points=coords, values=fn_value)
my_grid_interp(np.array([0, 0, 0]))  # Runs without error
my_rbf_interp = interp.Rbf(*[x.flatten() for x in coords_mesh], d=fn_value.flatten())
## Error: numpy.linalg.linalg.LinAlgError: singular matrix -- why?

What am I missing? The example above uses the function F(x, y, z) = x^2 + y*z. I'd like to use Rbf to approximate that function. As far as I can tell there are no duplicate coordinates: compare unique_coords_array to coords_array.

Upvotes: 5

Views: 2412

Answers (1)

armatita
armatita

Reputation: 13465

I believe the problem is your input:

    my_rbf_interp = interp.Rbf(*[x.flatten() for x in coords_mesh],d=fn_value.flatten())

Should you change to:

    x,y,z = [x.flatten() for x in coords_mesh]
    my_rbf_interp = interp.Rbf(x,y,z,fn_value.flatten())

And it should work. I think your original formulation is repeating lines in the matrix that goes for solve and thus having a very similar problem to duplicates (i.e. Singular Matrix).

Also if you would do:

    d = fn_value.flatten()
    my_rbf_interp = interp.Rbf(*(x,y,z,d))

It should work also.

Upvotes: 1

Related Questions