python - 这可能是 Anaconda 2.5.0 环境中 scipy.sparse.linalg.eigs 函数中的错误吗？

Question

我在 Linux Anaconda 2.5.0 环境中的 Jupyter notebook 中运行 scipy.sparse.linalg.eigs 函数，但每次运行时特征值似乎都在变化。当我在 Mac osX 环境中运行完全相同的笔记本时，无论我运行笔记本多少次，它都只为我提供一个特征值的解决方案，情况应该如此。这可能是 Linux 系统上 scipy 包的错误吗？我正在运行的笔记本如下:

import numpy as np
from scipy import fftpack as fft
from scipy.sparse import lil_matrix, csc_matrix
from scipy.sparse.linalg import eigs, inv

nz = 100
zin = np.arange(nz+1, dtype=np.float64)/nz
N2 = np.full(nz, 1.)
f0 = 1.
beta = 0.
Nx = int(1e2)
Ny = int(1e2)
dx = 1e-1
dy = 1e-1
vbar = np.zeros(nz+1)
ubar = zin
etax = np.zeros(2)
etay = np.zeros(2)
k = fft.fftfreq(Nx, dx)[:Nx/2]
l = fft.fftfreq(Ny, dy)[:Ny/2]

zc = np.hstack(0.5*(zin[1:] + zin[:-1]))
dzc = np.hstack(np.diff(zc))
depth = 1.
dztop = zc[0]
dzbot = depth - zc[-1]
zf = zin
dzf = np.diff(zf)

j = 0

for i in range(len(k)):

    L = lil_matrix((nz+1, nz+1), dtype=np.float64)
    G = L.copy()

    ################
    # n = 0 (surface)
    ################
    R = k[i] * .5*(ubar[0]+ubar[1]) + l[j] * .5*(vbar[0]+vbar[1]) 
    D = dzf[0]**-1
    S = .5 * ( k[i] * ( (ubar[0]-ubar[1])*D - N2[0]/f0 * etay[0] )
                    + l[j] * ( (vbar[0]-vbar[1])*D + N2[0]/f0 * etax[0] ) )

    L[0, 0] = R * D - S
    L[0, 1] = R * (-D) - S
    G[0, 0] = D
    G[0, 1] = - D

    ################
    # n = nz (bottom)
    ################
    R = k[i] * .5*(ubar[nz-1]+ubar[nz]) + l[j] * .5*(vbar[nz-1]+vbar[nz])
    D = dzf[nz-1]**-1
    S = .5 * ( k[i] * ( (ubar[nz-1]-ubar[nz])*D - N2[nz-1]/f0 * etay[1] )
                    + l[j] * ( (vbar[nz-1]-vbar[nz])*D + N2[nz-1]/f0 * etax[1] ) )

    L[nz, nz-1] = R * D - S
    L[nz, nz] = R * (-D) - S
    G[nz, nz-1] = D
    G[nz, nz] = - D

    ################
    # 0 < n < nz (interior)
    ################
    for n in range(1,nz):

        R = k[i] * ubar[n] + l[j] * vbar[n]
        K2 = k[i]**2 + l[j]**2
        bf = f0**2 * dzc[n-1]**-1
        b_1 = N2[n-1] * dzf[n-1]
        b = N2[n] * dzf[n]
        B_1 = bf * b_1**-1
        B = - (bf * ( b**-1 + b_1**-1 ) + K2)
        Bt1 = bf * b**-1

        N2Z = (N2[n]*dzf[n])**-1
        N2Z_1 = (N2[n-1]*dzf[n-1])**-1
        P = ( k[i] * ( beta - bf * ( ubar[n+1] * N2Z
                   - (N2Z + N2Z_1) * ubar[n]
                   + N2Z_1 * ubar[n-1] ) )
                   - l[j] * bf * ( vbar[n+1] * N2Z 
                   - (N2Z + N2Z_1) * vbar[n]
                   + N2Z_1 * vbar[n-1] )
            )

        L[n, n-1] = R * B_1
        L[n, n] = R * B + P
        L[n, n+1] = R * Bt1
        G[n, n-1] = B_1
        G[n, n] = B
        G[n, n+1] = Bt1

    val, func = eigs( csc_matrix(inv(csc_matrix(G)).dot(csc_matrix(L))), 
                                    k=2, which='LI', ncv=100, maxiter=1000 )  # default returns 6 eigenvectors

    if i == 0 and j == 0:
        omega = np.zeros( (len(val), len(k)), dtype=complex )
        psi = np.zeros( (nz+1, len(val), len(k)), dtype=complex )
    omega[:, i] = val
    psi[:, :, i] = func  # Each column is the eigenfunction

omega_imag = omega.imag.max(axis=1)
p = np.argsort(omega_imag)[::-1]
omega = omega[p]
psi = psi[:, p]

print(omega_imag[p])

我第一次在 Linux 环境中运行这段代码时，特征值(即 omega_imag)是

[ 0.30862832  0.28027858]

但是第二次给了

[ 0.3097162   0.26593829]

我运行的环境是:

# packages in environment at /home/takaya/.conda/envs/oceanmodes:
#
backports                 1.0                      py27_0    defaults
backports-abc             0.4                       <pip>
backports.shutil-get-terminal-size 1.0.0                     <pip>
backports.ssl-match-hostname 3.4.0.2                   <pip>
backports_abc             0.4                      py27_0    defaults
cairo                     1.12.18                       6    defaults
configparser              3.5.0b2                  py27_1    defaults
cycler                    0.10.0                   py27_0    defaults
decorator                 4.0.9                    py27_0    defaults
entrypoints               0.2                      py27_1    defaults
fontconfig                2.11.1                        5    defaults
freetype                  2.5.5                         0    defaults
functools32               3.2.3.2                  py27_0    defaults
get_terminal_size         1.0.0                    py27_0    defaults
ipykernel                 4.3.1                    py27_0    defaults
ipython                   4.2.0                    py27_0    defaults
ipython-genutils          0.1.0                     <pip>
ipython_genutils          0.1.0                    py27_0    defaults
ipywidgets                4.1.1                    py27_0    defaults
jinja2                    2.8                      py27_0    defaults
jsonschema                2.5.1                    py27_0    defaults
jupyter                   1.0.0                    py27_2    defaults
jupyter-client            4.2.2                     <pip>
jupyter-console           4.1.1                     <pip>
jupyter-core              4.1.0                     <pip>
jupyter_client            4.2.2                    py27_0    defaults
jupyter_console           4.1.1                    py27_0    defaults
jupyter_core              4.1.0                    py27_0    defaults
libgcc                    5.2.0                         0    defaults
libgfortran               3.0.0                         1    defaults
libpng                    1.6.17                        0    defaults
libsodium                 1.0.3                         0    defaults
libxml2                   2.9.2                         0    defaults
markupsafe                0.23                     py27_0    <unknown>
matplotlib                1.5.1               np111py27_0    defaults
mistune                   0.7.2                    py27_0    defaults
mkl                       11.3.1                        0    defaults
nbconvert                 4.2.0                    py27_0    defaults
nbformat                  4.0.1                    py27_0    defaults
notebook                  4.2.0                    py27_0    defaults
numpy                     1.11.0                   py27_0    defaults
openssl                   1.0.2h                        0    defaults
path.py                   8.2.1                    py27_0    defaults
pexpect                   4.0.1                    py27_0    defaults
pickleshare               0.5                      py27_0    defaults
pip                       8.1.1                    py27_1    defaults
pixman                    0.32.6                        0    defaults
ptyprocess                0.5                      py27_0    defaults
pycairo                   1.10.0                   py27_0    defaults
pygments                  2.1.3                    py27_0    defaults
pyparsing                 2.1.1                    py27_0    defaults
pyqt                      4.11.4                   py27_1    defaults
python                    2.7.11                        0    defaults
python-dateutil           2.5.2                    py27_0    defaults
pytz                      2016.3                   py27_0    defaults
pyzmq                     15.2.0                   py27_0    defaults
qt                        4.8.7                         1    defaults
qtconsole                 4.2.1                    py27_0    defaults
readline                  6.2                           2    <unknown>
scipy                     0.17.0              np111py27_3    defaults
setuptools                20.7.0                   py27_0    defaults
simplegeneric             0.8.1                    py27_0    defaults
singledispatch            3.4.0.3                  py27_0    defaults
sip                       4.16.9                   py27_0    defaults
six                       1.10.0                   py27_0    defaults
sqlite                    3.9.2                         0    defaults
ssl_match_hostname        3.4.0.2                  py27_1    defaults
terminado                 0.5                      py27_1    defaults
tk                        8.5.18                        0    http://repo.continuum.io/pkgs/free/linux-64/tk-8.5.18-0.tar.bz2
tornado                   4.3                      py27_0    defaults
traitlets                 4.2.1                    py27_0    defaults
wheel                     0.29.0                   py27_0    defaults
zeromq                    4.1.3                         0    defaults
zlib                      1.2.8                         0    http://repo.continuum.io/pkgs/free/linux-64/zlib-1.2.8-0.tar.bz2

Answer 1

ARPACK中的起始向量，如果没有指定，是随机的。您可以尝试指定 v0 参数以使用固定参数。

此外，如果计算对初始值敏感(病态​​)，则可能会在不同的运行中得到不同的结果，因为浮点舍入误差会受到可能变化的外部因素(内存)的影响对齐、BLAS 中的线程等)。

换句话说，如今要确保确定性结果需要的不仅仅是编写并运行程序——您可能需要使用额外的 C/Fortran 编译器选项来编译所有内容(Numpy、scipy 等)，以产生更严格的确定性 float -点运算机器码。您可能还需要在优化的线性代数库(MKL、Openblas)中控制线程等的使用。最简单的选择通常是确保计算不是病态的，这样小的扰动就无关紧要了——但在 ARPACK 的情况下，这可能不在你的直接控制之下，除非它可能以这样的变化无关紧要(或者您根本不计算特征值)的方式构建您尝试进行的计算。