SVD did not converge

[Ml0749276]

I got this while I tried to OrdinaryKriging in a loop.

Here is part of code

        for obse in ObseSites.items:
            for j in range(len(obse.freqs)):
                if math.isclose(obse.freqs[j], freq):
                    x.append(obse.x/1000.)
                    y.append(obse.y/1000.)
                    rxy.append(math.log10(obse.rxy[j].re))
                    pxy.append(DegtoFirstQuadrant(obse.pxy[j].re))
                    ryx.append(math.log10(obse.ryx[j].re))
                    pyx.append(DegtoFirstQuadrant(obse.pyx[j].re))
                    break

        for j in range(len(component)):
            ok = OrdinaryKriging(y, x, component[j], variogram_model='exponential', verbose=True)
            z, ss = ok.execute('grid', gridy, gridx)

The error output is as below

Adjusting data for anisotropy...
Initializing variogram model...

[MuellerSeb]

Hey there!

I got this while I tried to OrdinaryKriging in a loop.

Where did it say something about SVD as you have written in the title?

The error output is as below

Adjusting data for anisotropy...
Initializing variogram model...

Is there really an error? Or is it hanging at this point? If it is hanging, I think it could be at this line:

PyKrige/pykrige/core.py

Line 455 in b33e5d8

d = pdist(X, metric="euclidean")

Or did you miss to copy the full output? How big is the size of x and y in you case?

Cheers, Sebastian

[Ml0749276]

Hi Sebastian,

Here is the full script

# filecoding=utf-8

# 用于对测点数据进行插值，得到平面值

import math
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
from pykrige.ok import OrdinaryKriging
import readZdatafile
import sitesclass

D2R = math.pi/180
R2D = 180/math.pi


def DegtoFirstQuadrant(angle):
    result = math.atan(math.tan(angle*D2R))*R2D
    if result < 0:
        result = result+180
    return result


Tsite = sitesclass.Tsite
Tsites = sitesclass.Tsites
ObseSites = Tsites()


obsesitelist = []
obsesitelist.append('0m0d5%.ori_dat')
obsesitelist.append('0m85d5%.ori_dat')


# 待绘制频率
freqlist = [0.25, 0.05, 0.025, 0.0125, 0.00625]


for i in range(len(obsesitelist)):
    obsedir = obsesitelist[i]
    readZdatafile.readfiles(obsedir, ObseSites)

    # 设置网格
    xstep = (ObseSites.xmax-ObseSites.xmin+20)/(len(ObseSites.items)-1)
    ystep = (ObseSites.ymax-ObseSites.ymin+20)/(len(ObseSites.items)-1)
    gridx = np.arange(ObseSites.xmin/1000.-10, ObseSites.xmax/1000.+10, min(xstep, ystep)/1000.)
    gridy = np.arange(ObseSites.ymin/1000.-10, ObseSites.ymax/1000.+10, min(xstep, ystep)/1000.)

    fig, ax = plt.subplots(4, 5, dpi=180, sharex='all', sharey='all', figsize=(4, 5))

    # 循环频率，绘制观测数据
    for index in range(len(freqlist)):
        freq = freqlist[index]
        x = np.zeros(len(ObseSites.items))
        y = np.zeros(len(ObseSites.items))
        rxy = np.zeros(len(ObseSites.items))
        pxy = np.zeros(len(ObseSites.items))
        ryx = np.zeros(len(ObseSites.items))
        pyx = np.zeros(len(ObseSites.items))
        component = [rxy, ryx, pxy, pyx]

        for obse in ObseSites.items:
            iobse = ObseSites.items.index(obse)
            for j in range(len(obse.freqs)):
                if math.isclose(obse.freqs[j], freq):
                    x[iobse] = obse.x/1000.
                    y[iobse] = obse.y/1000.
                    rxy[iobse] = math.log10(obse.rxy[j].re)
                    pxy[iobse] = DegtoFirstQuadrant(obse.pxy[j].re)
                    ryx[iobse] = math.log10(obse.ryx[j].re)
                    pyx[iobse] = DegtoFirstQuadrant(obse.pyx[j].re)
                    break

        for j in range(len(component)):
            if j < 2:
                vmin = 1
                vmax = 5
            else:
                vmin = 0
                vmax = 100

            ok = OrdinaryKriging(y, x, component[j], variogram_model='exponential', verbose=True)
            z, ss = ok.execute('grid', gridy, gridx)

            ax[j, index].pcolormesh(gridy, gridx, z, cmap='jet_r',
                                    vmin=vmin, vmax=vmax, linewidth=0)
            ax[j, index].scatter(y, x, c='black', s=5)
            ax[j, index].set_aspect(1)
            ax[j, index].tick_params(labelsize=5)
            # 调整坐标轴
            ax[j, 0].set_ylabel('X(km)', fontsize=5)

        # 设置图头
        ax[0, index].set_title(str(freq)+'Hz', fontsize=5, pad=0.6)
        # 调整坐标轴
        ax[3, index].set_xlabel('Y(km)', fontsize=5)

    # 调整布局
    fig.subplots_adjust(hspace=0.1, wspace=0.05, right=0.85,
                        left=0.11, top=0.95, bottom=0.1)
    # 设置右侧分量名
    fig.text(0.85, 0.85, 'Rxy', fontsize=5, va='center', rotation=90)
    fig.text(0.85, 0.63, 'Ryx', fontsize=5, va='center', rotation=90)
    fig.text(0.85, 0.41, 'Pxy', fontsize=5, va='center', rotation=90)
    fig.text(0.85, 0.20, 'Pyx', fontsize=5, va='center', rotation=90)
    # 设置色标
    # 视电阻率
    norm = mpl.colors.Normalize(vmin=1, vmax=5)
    m = cm.ScalarMappable(norm=norm, cmap=cm.jet_r)
    pos = fig.add_axes([0.9, 0.55, 0.02, 0.35])  # 位置[左、下、宽度、高度]
    cb = plt.colorbar(m, cax=pos)
    cb.ax.tick_params(labelsize=5)
    cb.ax.set_title('log10(rho)', fontsize=5, loc='left')
    # 相位
    norm = mpl.colors.Normalize(vmin=0, vmax=100)
    m = cm.ScalarMappable(norm=norm, cmap=cm.jet_r)
    pos = fig.add_axes([0.9, 0.11, 0.02, 0.35])  # 位置[左、下、宽度、高度]
    cb = plt.colorbar(m, cax=pos)
    cb.ax.tick_params(labelsize=5)
    cb.ax.set_title('degree', fontsize=5, loc='left')

    plt.show()
    plt.close()

I write this script to genarate figures, the size of 'x' and 'y' is 33, and the size of the mesh is [38,87]. During the whole workflow, I read two data files and generate two figures for each one, after the first figure generated successfully, the second one will jump out the error 'SVD did not converge' as attached picuture. In the beginning, I thought it might be the problem of my datafile, so I switched the order of the two datafile and read the 2nd first, then it worked and generated figure for the 2nd data file, but the same error occured while processing the 1st data file. So I got confused. Can you help me solving this? Thank you!

Attched are the error message captured and two figures generated. PyKrige is a convenient tool for the interpolation and plotting!

The error message

The 1st data file's figure

The 2nd data file's figure

[Ml0749276]

And another question, this subfigure seems not interpolate correctly, do you have some ideas about this? I can send you all the data if you want, this situation not only appeared here, but many other places.

[MuellerSeb]

Also in you case it could help to use the pseudo-inverse version from this branch:
https://github.com/GeoStat-Framework/PyKrige/tree/pseudo_inv

Then you can alter your call:

ok = OrdinaryKriging(y, x, component[j], variogram_model='exponential', verbose=True, pseudo_inv=True)

Maybe this helps. This will be in the next version.

[Ml0749276]

Hi Sebastian,

Thanks for your reply, I have tried as you said, it makes no difference.

[MuellerSeb]

Could you check your data for nans and infs? That is often related to a svd problem

[Ml0749276]

Hi Sebastian,

I am on a business trip these days, I check the data and no nan or inf is found.