2014-12-20 11 views
12

Tôi có một tập lệnh không sử dụng phép ngẫu nhiên nào mang lại cho tôi các câu trả lời khác nhau khi tôi chạy nó. Tôi hy vọng câu trả lời sẽ giống nhau, mỗi khi tôi chạy kịch bản. Vấn đề dường như chỉ xảy ra đối với một số dữ liệu đầu vào (không có điều kiện).Kịch bản lệnh python xác định hoạt động theo cách không xác định

Đoạn mã xuất phát từ thuật toán để tính toán loại bộ điều khiển cụ thể cho hệ thống tuyến tính và chủ yếu bao gồm thực hiện đại số tuyến tính (ma trận nghịch đảo, phương trình Riccati, giá trị riêng).

Rõ ràng, đây là một lo lắng lớn đối với tôi, vì bây giờ tôi không thể tin tưởng mã của tôi để cho tôi kết quả đúng. Tôi biết kết quả có thể sai đối với dữ liệu có điều kiện kém, nhưng tôi mong đợi sai một cách nhất quán. Tại sao câu trả lời không phải lúc nào cũng giống nhau trên máy tính Windows của tôi? Tại sao máy tính Windows & không cho kết quả tương tự?

Tôi đang sử dụng Python 2.7.9 (default, Dec 10 2014, 12:24:55) [MSC v.1500 32 bit (Intel)] on win 32, với phiên bản Numpy 1.8.2 và Scipy 0.14.0. (Windows 8, 64bit).

Mã bên dưới. Tôi cũng đã thử chạy mã trên hai máy Linux, và có kịch bản luôn luôn đưa ra câu trả lời tương tự (nhưng các máy đã đưa ra các câu trả lời khác nhau). Một người đang chạy Python 2.7.8, với Numpy 1.8.2 và Scipy 0.14.0. Thứ hai là chạy Python 2.7.3 với Numpy 1.6.1 và Scipy 0.12.0.

Tôi giải phương trình Riccati ba lần, rồi in câu trả lời. Tôi hy vọng cùng một câu trả lời mỗi lần, thay vào đó tôi nhận được chuỗi '1.75305103767e-09; 3.25501787302e-07; 3.25501787302e-07 '.

import numpy as np 
    import scipy.linalg 

    matrix = np.matrix 

    A = matrix([[ 0.00000000e+00, 2.96156260e+01, 0.00000000e+00, 
         -1.00000000e+00], 
        [ -2.96156260e+01, -6.77626358e-21, 1.00000000e+00, 
         -2.11758237e-22], 
        [ 0.00000000e+00, 0.00000000e+00, 2.06196064e+00, 
         5.59422224e+01], 
        [ 0.00000000e+00, 0.00000000e+00, 2.12407340e+01, 
         -2.06195974e+00]]) 
    B = matrix([[  0.  ,  0.  ,  0.  ], 
        [  0.  ,  0.  ,  0.  ], 
        [ -342.35401351, -14204.86532216,  31.22469724], 
        [ 1390.44997337, 342.33745324, -126.81720597]]) 
    Q = matrix([[ 5.00000001, 0.  , 0.  , 0.  ], 
        [ 0.  , 5.00000001, 0.  , 0.  ], 
        [ 0.  , 0.  , 0.  , 0.  ], 
        [ 0.  , 0.  , 0.  , 0.  ]]) 
    R = matrix([[ -3.75632852e+04, -0.00000000e+00, 0.00000000e+00], 
        [ -0.00000000e+00, -3.75632852e+04, 0.00000000e+00], 
        [ 0.00000000e+00, 0.00000000e+00, 4.00000000e+00]]) 

    counter = 0 
    while counter < 3: 
      counter +=1 

      X = scipy.linalg.solve_continuous_are(A, B, Q, R) 
      print(-3449.15531628 - X[0,0]) 

cấu hình NumPy của tôi là như sau print np.show_config()

 
lapack_opt_info: 
    libraries = ['mkl_blas95', 'mkl_lapack95', 'mkl_intel_c', 'mkl_intel_thread', 'mkl_core', 'libiomp5md', 'mkl_blas95', 'mkl_lapack95', 'mkl_intel_c', 'mkl_intel_thread', 'mkl_core', 'libiomp5md'] 
    library_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/lib/ia32', 'C:/Program Files (x86)/Intel/Composer XE 2013 SP1/compiler/lib/ia32'] 
    define_macros = [('SCIPY_MKL_H', None)] 
    include_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/include'] 
blas_opt_info: 
    libraries = ['mkl_blas95', 'mkl_lapack95', 'mkl_intel_c', 'mkl_intel_thread', 'mkl_core', 'libiomp5md'] 
    library_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/lib/ia32', 'C:/Program Files (x86)/Intel/Composer XE 2013 SP1/compiler/lib/ia32'] 
    define_macros = [('SCIPY_MKL_H', None)] 
    include_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/include'] 
openblas_info: 
    NOT AVAILABLE 
lapack_mkl_info: 
    libraries = ['mkl_blas95', 'mkl_lapack95', 'mkl_intel_c', 'mkl_intel_thread', 'mkl_core', 'libiomp5md', 'mkl_blas95', 'mkl_lapack95', 'mkl_intel_c', 'mkl_intel_thread', 'mkl_core', 'libiomp5md'] 
    library_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/lib/ia32', 'C:/Program Files (x86)/Intel/Composer XE 2013 SP1/compiler/lib/ia32'] 
    define_macros = [('SCIPY_MKL_H', None)] 
    include_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/include'] 
blas_mkl_info: 
    libraries = ['mkl_blas95', 'mkl_lapack95', 'mkl_intel_c', 'mkl_intel_thread', 'mkl_core', 'libiomp5md'] 
    library_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/lib/ia32', 'C:/Program Files (x86)/Intel/Composer XE 2013 SP1/compiler/lib/ia32'] 
    define_macros = [('SCIPY_MKL_H', None)] 
    include_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/include'] 
mkl_info: 
    libraries = ['mkl_blas95', 'mkl_lapack95', 'mkl_intel_c', 'mkl_intel_thread', 'mkl_core', 'libiomp5md'] 
    library_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/lib/ia32', 'C:/Program Files (x86)/Intel/Composer XE 2013 SP1/compiler/lib/ia32'] 
    define_macros = [('SCIPY_MKL_H', None)] 
    include_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/include'] 
None 

(chỉnh sửa để cắt câu hỏi xuống)

+0

Bạn có gieo hạt máy phát điện không? http://docs.scipy.org/doc/numpy/reference/generated/numpy.random.seed.html – NPE

+0

Tôi đã thêm 'np.random.seed (0)' làm dòng mã đầu tiên (ngay sau đầu vào), và nó không có sự khác biệt. Sẽ cập nhật câu hỏi. – mwmwm

+0

Xin lỗi, tôi bằng cách nào đó đã hiểu sai câu hỏi của bạn vì nói rằng mã * không * sử dụng ngẫu nhiên. Rất tiếc. – NPE

Trả lời

7

Nói chung, thư viện linalg trên Windows cung cấp cho câu trả lời khác nhau trên chạy khác nhau ở mức độ chính xác máy . Tôi chưa bao giờ nghe giải thích tại sao điều này chỉ xảy ra hoặc chủ yếu trên Windows.

Nếu ma trận của bạn bị ảnh hưởng, thì cuộc gọi sẽ phần lớn là nhiễu số. Trên Windows, tiếng ồn không phải lúc nào cũng giống nhau trong các lần chạy liên tiếp, trên các hệ điều hành khác, nhiễu có thể luôn giống nhau nhưng có thể khác nhau tùy thuộc vào chi tiết của thư viện đại số tuyến tính, tùy chọn luồng, sử dụng bộ nhớ cache và vân vân.

Tôi đã xem và đăng lên danh sách gửi thư scipy một vài ví dụ về điều này trên Windows, tôi đã sử dụng hầu hết các tệp nhị phân 32 bit chính thức với ATLAS BLAS/LAPACK.

Giải pháp duy nhất là làm cho kết quả tính toán của bạn không phụ thuộc quá nhiều vào vấn đề độ chính xác điểm nổi và nhiễu số, ví dụ như chuẩn hóa ma trận nghịch đảo, sử dụng nghịch đảo tổng quát, pinv, reparameterize hoặc tương tự.

+3

Một phần của nhiễu số xuất phát từ liên kết dữ liệu trong bộ nhớ, tức là liệu nó có xảy ra trực tiếp ở vị trí phù hợp với lệnh sse vv hay không --- hai trường hợp hoạt động theo thứ tự khác nhau, do đó lỗi làm tròn điểm nổi là khác nhau. win32 tiêu chuẩn bộ nhớ cấp phát afaik có xu hướng trở lại khối bộ nhớ được trên trung bình ít liên kết hơn phân bổ Linux, do đó, có jitter hơn ở đó. Nếu kết quả của bạn phụ thuộc đáng kể vào cách mà lỗi làm tròn FP đi, thì việc xây dựng vấn đề của bạn không ổn định về mặt số lượng. –

2

pv được ghi chú trong comments to user333700's answer, công thức trước đây của trình giải quyết Riccati là, mặc dù về lý thuyết chính xác, không ổn định về số lượng. Vấn đề này được khắc phục trên phiên bản phát triển của SciPy và những người giải quyết hỗ trợ phương trình Riccati tổng quát.

Trình giải quyết Riccati được cập nhật và kết quả giải quyết sẽ có sẵn từ phiên bản 0.19 trở đi. Bạn có thể kiểm tra số development branch docs here.

Nếu, bằng cách sử dụng ví dụ được đưa ra trong câu hỏi tôi thay thế các vòng cuối cùng với

for _ in range(5): 
    x = scipy.linalg.solve_continuous_are(A, B, Q, R) 
    Res = [email protected] + [email protected] + q - [email protected]@ np.linalg.solve(r,b.T)@ x 
    print(Res) 

tôi nhận được (cửa sổ 10, py3.5.2)

[[ 2.32314924e-05 -2.55086270e-05 -7.66709854e-06 -9.01878229e-06] 
[ -2.62447211e-05 2.61182140e-05 8.27328768e-06 1.00345896e-05] 
[ -7.92257197e-06 8.57094892e-06 2.50908488e-06 3.05714639e-06] 
[ -9.51046241e-06 9.80847472e-06 3.13103374e-06 3.60747799e-06]] 
[[ 2.32314924e-05 -2.55086270e-05 -7.66709854e-06 -9.01878229e-06] 
[ -2.62447211e-05 2.61182140e-05 8.27328768e-06 1.00345896e-05] 
[ -7.92257197e-06 8.57094892e-06 2.50908488e-06 3.05714639e-06] 
[ -9.51046241e-06 9.80847472e-06 3.13103374e-06 3.60747799e-06]] 
[[ 2.32314924e-05 -2.55086270e-05 -7.66709854e-06 -9.01878229e-06] 
[ -2.62447211e-05 2.61182140e-05 8.27328768e-06 1.00345896e-05] 
[ -7.92257197e-06 8.57094892e-06 2.50908488e-06 3.05714639e-06] 
[ -9.51046241e-06 9.80847472e-06 3.13103374e-06 3.60747799e-06]] 
[[ 2.32314924e-05 -2.55086270e-05 -7.66709854e-06 -9.01878229e-06] 
[ -2.62447211e-05 2.61182140e-05 8.27328768e-06 1.00345896e-05] 
[ -7.92257197e-06 8.57094892e-06 2.50908488e-06 3.05714639e-06] 
[ -9.51046241e-06 9.80847472e-06 3.13103374e-06 3.60747799e-06]] 
[[ 2.32314924e-05 -2.55086270e-05 -7.66709854e-06 -9.01878229e-06] 
[ -2.62447211e-05 2.61182140e-05 8.27328768e-06 1.00345896e-05] 
[ -7.92257197e-06 8.57094892e-06 2.50908488e-06 3.05714639e-06] 
[ -9.51046241e-06 9.80847472e-06 3.13103374e-06 3.60747799e-06]] 

Để tham khảo, giải pháp trở lại là

array([[-3449.15531305, 4097.1707738 , 1142.52971904, 1566.51333847], 
     [ 4097.1707738 , -4863.70533241, -1356.66524959, -1860.15980716], 
     [ 1142.52971904, -1356.66524959, -378.32586814, -518.71965102], 
     [ 1566.51333847, -1860.15980716, -518.71965102, -711.21062988]]) 
Các vấn đề liên quan