在 Python 的输出中“nan”而不是矩阵 elemets?

"nan" instead of matrix elemets in output in Python?

提问人:Anushka Tilekar 提问时间:11/5/2023 最后编辑:Anushka Tilekar 更新时间:11/6/2023 访问量:58

问:

我正在尝试计算这个稀疏矩阵问题 Au=f 的近似解,但我在矩阵的输出中得到“nan”而不是矩阵元素,尽管输入矩阵都包含 np.float64。 我不明白我哪里出了问题。

这是我的代码:-

import numpy as np
import scipy
from scipy.sparse import coo_matrix
from scipy.sparse import random

#Defining a function that takes N as an input

def our_func (N) :

  #Taking the given values of parameters h and k
  h =  1/N
  k = (29*(np.pi))/2
   
  #Defining initial input matrix with domensions N+1
  I = np.empty([N+1,N+1],  dtype=np.float64) 
  #[REF :- https://www.geeksforgeeks.org/how-to-create-an-empty-matrix-with-numpy-in-python/]
  print("datatype of I is :-")
  print(I.dtype)

  #Generating the matrix A in coo sparse format
  A_coo = coo_matrix(I).reshape((N+1,N+1))
  print("datatype of A_coo is :-")
  print(A_coo.dtype)

  #Converting the matrix A from coo sparse format to CSR sparse format 
  A = A_coo.tocsr().reshape((N+1,N+1))
  print("datatype of A is :-")
  print(A.dtype)
 
  # Creating an empty storage array for the right side of our equation Au=f
  f = np.empty((N+1), dtype=np.float64)

  #The rows of matrix A and f are defined as :-
  A[0,0] = 1    #for i==j==0
  A[N,N] = 1    #for i==j==N
  f[N] = 1      #for i==N
  f[0] = 0      #for i==0
  for i in range (int(N)) :
    for j in range (int(N)) :
        A[i,1] = 2 - ((h**2)*(k**2)) 
        A[i,i+1] = -1
        A[i,i-1] = -1
        A[i,j] = 0      #for i!==j
  
  #Our function Returns matrix A and f, with matrix A in sparse storage format
  return A, f

# Defining a new function using the function scipy.sparse.linalg.spsolve 
# to solve the sparse matrix-vector problem: Au=f

def our_solution(N):
  A,f = our_func(N)
  our_sparse_matrix = scipy.sparse.linalg.spsolve(A,f)
  return our_sparse_matrix

# Using this defined function to compute the approximate solution for our 
# problem, for N=10, N=100, and N=1000.
approx_sol_1 = our_solution(10)
print("approx_sol_1 is :- ")
print(approx_sol_1)

python numpy scipy nan 稀疏矩阵

评论

0赞 hpaulj 11/5/2023
我不知道这是否重要,但这是一个失败。nan

答:

0赞 hpaulj 11/5/2023 #1

为较小的函数运行函数:N

In [4]: A,f = our_func(10)
datatype of I is :-
float64
datatype of A_coo is :-
float64
datatype of A is :-
float64
C:\Users\14256\miniconda3\lib\site-packages\scipy\sparse\_index.py:103: SparseEfficiencyWarning: Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient.
  self._set_intXint(row, col, x.flat[0])

它警告我们,单独设置稀疏数组的元素是没有效率的。csr

In [5]: A
Out[5]: 
<11x11 sparse matrix of type '<class 'numpy.float64'>'
    with 103 stored elements in Compressed Sparse Row format>

注意的是,矩阵不是特别稀疏,103 个值中有 121 个值是非零的。

In [6]: A.A
Out[6]: 
array([[  0.        ,  -1.        ,   0.        ,   0.        ,
          0.        ,   0.        ,   0.        ,   0.        ,
          0.        ,   0.        ,  -1.        ],
       [ -1.        , -18.75084325,  -1.        ,   0.        ,
          0.        ,   0.        ,   0.        ,   0.        ,
          0.        ,   0.        ,   0.        ],...

它可以被修剪,摆脱许多 0。

In [10]: A.eliminate_zeros()
In [11]: A
Out[11]: 
<11x11 sparse matrix of type '<class 'numpy.float64'>'
    with 28 stored elements in Compressed Sparse Row format>

In [17]: def my_func (N) :
    ...: 
    ...:   #Taking the given values of parameters h and k
    ...:   h =  1/N
    ...:   k = (29*(np.pi))/2
    ...: 
    ...:   #Defining initial input matrix with domensions N+1
    ...:   A = np.zeros([N+1,N+1],  dtype=np.float64)
    ...: 
    ...:   # Creating an empty storage array for the right side of our equation Au=f
    ...:   f = np.zeros((N+1), dtype=np.float64)
    ...: 
    ...:   #The rows of matrix A and f are defined as :-
    ...:   A[0,0] = 1    #for i==j==0
    ...:   A[N,N] = 1    #for i==j==N
    ...:   f[N] = 1      #for i==N
    ...:   #f[0] = 0      #for i==0
    ...:   for i in range (int(N)) :
    ...:     for j in range (int(N)) :
    ...:         A[i,1] = 2 - ((h**2)*(k**2)) 
    ...:         A[i,i+1] = -1
    ...:         A[i,i-1] = -1
    ...:         A[i,j] = 0      #for i!==j
    ...: 
    ...:   #Our function Returns matrix A and f, with matrix A in sparse storage format
    ...:   return A, f
    ...:   

In [18]: A,f = my_func(10)

In [19]: A.shape
Out[19]: (11, 11)

这可以在以下情况下转换为稀疏:

In [21]: M = coo_matrix(A).tocsr()
In [22]: M
Out[22]: 
<11x11 sparse matrix of type '<class 'numpy.float64'>'
    with 30 stored elements in Compressed Sparse Row format>

除非太大了,否则你会遇到内存错误,我认为这种方法会更快。N

这并不能解决您的价值观问题。我对你的问题了解不够,不知道是否有合理的值。nanA

编辑

我想知道是不是这样

A[i,j] = 0      #for i!==j

太宽泛了。如果我在 中省略它,推断它已经是 0,则有一个不同的值 (-1)。my_funcA[8,9]

In [51]: A,f = my_func(10)

In [52]: scipy.sparse.linalg.spsolve(A,f)
C:\Users\14256\miniconda3\lib\site-packages\scipy\sparse\linalg\_dsolve\linsolve.py:214: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format
  warn('spsolve requires A be CSC or CSR matrix format',
Out[52]: 
array([ 0.8987385, -0.1012615,  1.       ,  0.1012615,  0.8987385,
        1.797477 ,  1.       ,  0.1012615,  0.8987385,  1.797477 ,
        1.       ])

使用自:A1our_func

In [53]: scipy.sparse.linalg.spsolve(A1,f)
C:\Users\14256\miniconda3\lib\site-packages\scipy\sparse\linalg\_dsolve\linsolve.py:276: MatrixRankWarning: Matrix is exactly singular
  warn("Matrix is exactly singular", MatrixRankWarning)
Out[53]: array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan])

更正稀疏警告会给出相同的结果:

In [54]: scipy.sparse.linalg.spsolve(coo_matrix(A).tocsr(),f)
Out[54]: 
array([ 0.8987385, -0.1012615,  1.       ,  0.1012615,  0.8987385,
        1.797477 ,  1.       ,  0.1012615,  0.8987385,  1.797477 ,
        1.       ])

我不知道这是否正确,至少它不是单一的。看起来应该有 3 个对角线,“-1,?, -1”。AA

0赞 Matt Haberland 11/6/2023 #2

我不明白我哪里出了问题。

矩阵是奇异的,因此您尝试求解的线性系统没有唯一的解。您可以看到这一点,因为倒数第二列的所有元素都为 和 为零(或基本上为零)。例如:AN=10N=100

A, f = our_func(10)
A.todense()[:, -2].T
# matrix([[0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
#          0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
#          0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
#          0.00000000e+000, 4.69807001e-294]])
# 4.69807001e-294 is essentially 0

当遇到(数值)奇异矩阵时,它会发出警告并返回 NaN 数组,而不是引发错误。sparse.linalg.spsolve

from scipy import sparse
A = sparse.csc_matrix((3, 3))
b = np.zeros(3)
sparse.linalg.spsolve(A, b)
# array([nan, nan, nan])

要解决此问题,请检查用于生成矩阵的代码,并确保其按预期工作。如果您确定这是您尝试求解的线性系统,则需要确定系统没有唯一解意味着什么。