Array changes unexpectedly when used with numpy.linalg functions

[walkerdray]

I have a jitclass where I want to save an eigenvector matrix. Later, I use that eigenvector matrix to solve some linear systems. However, when passed to numpy's linalg.solve function, the original eigenvector matrix is changed in place. Here's a minimal example:

import numpy as np

import numba as nb
from numba import njit, jitclass

from collections import OrderedDict

spec = OrderedDict(G=nb.float64[:,:])

@jitclass(spec)
class cls_jit(object):
    def __init__(self):
        pass
        
        
    def test0(self):
        x = np.ones(2)
        M = np.array([[1.,2.],[2.,1.]])
        _, G = np.linalg.eig(M)
        self.G = G
        save_G = []
        save_G.append(self.G.copy())
        for i in range(4):
            np.linalg.solve(self.G, x)
            save_G.append(self.G.copy())
        return save_G
    
    def test1(self):
        x = np.ones(2)
        M = np.array([[1.,2.],[2.,1.]])
        _, G = np.linalg.eig(M)
        save_G = []
        save_G.append(G.copy())
        for i in range(4):
            np.linalg.solve(G, x)
            save_G.append(G.copy())
        return save_G

@njit
def test2():
    x = np.ones(2)
    M = np.array([[1.,2.],[2.,1.]])
    _, G = np.linalg.eig(M)
    save_G = []
    save_G.append(G.copy())
    for i in range(4):
        np.linalg.solve(G, x)
        save_G.append(G.copy())
    return save_G

CJ = cls_jit()

print('running test0')
save_g = CJ.test0()
for g in save_g:
    print(g)
    print()
# output:
# running test0
# [[ 0.70710678 -0.70710678]
#  [ 0.70710678  0.70710678]]
# 
# [[ 0.70710678 -0.70710678]
#  [ 1.          1.41421356]]
# 
# [[ 1.          1.41421356]
#  [ 0.70710678 -1.70710678]]
# 
# [[ 1.          1.41421356]
#  [ 0.70710678 -2.70710678]]
# 
# [[ 1.          1.41421356]
#  [ 0.70710678 -3.70710678]]

print('running test1')
save_g = CJ.test1()
for g in save_g:
    print(g)
    print()
# output:
# running test1
# [[ 0.70710678 -0.70710678]
#  [ 0.70710678  0.70710678]]
# 
# [[ 0.70710678 -0.70710678]
#  [ 0.70710678  0.70710678]]
# 
# [[ 0.70710678 -0.70710678]
#  [ 0.70710678  0.70710678]]
# 
# [[ 0.70710678 -0.70710678]
#  [ 0.70710678  0.70710678]]
# 
# [[ 0.70710678 -0.70710678]
#  [ 0.70710678  0.70710678]]

print('running test2')
save_g = test2()
for g in save_g:
    print(g)
    print()
# output:
# running test2
# [[ 0.70710678 -0.70710678]
#  [ 0.70710678  0.70710678]]
# 
# [[ 0.70710678 -0.70710678]
#  [ 0.70710678  0.70710678]]
# 
# [[ 0.70710678 -0.70710678]
#  [ 0.70710678  0.70710678]]
# 
# [[ 0.70710678 -0.70710678]
#  [ 0.70710678  0.70710678]]
# 
# [[ 0.70710678 -0.70710678]
#  [ 0.70710678  0.70710678]]

I'm using python 3.6, numba version 0.45.1 and numpy version 1.16.4. In the above code, test1 and test2 work as expected. However, in test0, the value of G changes after every call to np.linalg.solve.

I'm guessing this has something to do with the fact that np.linalg.eig returns the eigenvectors as a view, and this somehow interacts with being assigned to the jitclass.

[sklam]

I can replicate this bug. This is definitely jitclass related because the problem doesn't occur if G is stored and retrieved from typed-list or typed-dict.

Here's a smaller reproducer:

import numpy as np

import numba as nb
from numba import njit, jitclass

from collections import OrderedDict

spec = OrderedDict(G=nb.float64[:,:])

@jitclass(spec)
class cls_jit(object):
    def __init__(self):
        pass


@njit
def test0(self, M, x):
    _, G = np.linalg.eig(M)
    self.G = G
    for i in range(4):
        A = self.G
        print('--', np.linalg.norm(A))  # norm of A should not change
        np.linalg.solve(A, x) 


CJ = cls_jit()

print('running test0')
x = np.ones(2)
M = np.array([[1.,2.],[2.,1.]])
test0(CJ, M, x)

[stuartarchibald]

Pretty sure this is a bug in linalg that I've seen before, not jitclass related, it's just that the type system makes it easy(-y/ier) to hit it. Reproducer that has no jitclass involved:

import numpy as np
from numba import njit, float64

@njit(locals={'G': float64[:,:]})
def foo(A):
    _, G = np.linalg.eig(A)
    for i in range(3):
        np.linalg.solve(G, G[:, 0])
        print('--', np.linalg.norm(G))

M = np.array([[1.,2.],[2.,1.]])
foo(M)

I think the problem is that if an array is determined to be 'A' order at typing time but at runtime it is discovered as 'F', then no copy is made:

numba/numba/targets/linalg.py

Lines 1721 to 1725 in 60d2bdd

	# a is destroyed on exit, copy it
	if a_F_layout:
	acpy = np.copy(a)
	else:
	acpy = np.asfortranarray(a)

[stuartarchibald]

xref #3368

[stuartarchibald]

Fixed by #5879. Closing.