آزمون Cholesky کوواریانس با NumPy؛ تقارن، نیمهٔ ماتریس و عامل مثلثی

یک فراخوانیِ موفقِ تجزیهٔ ماتریس، لزوماً اعتبارِ همهٔ داده‌های ورودی را ثابت نمی‌کند. در بازبینیِ ماژول کوواریانسِ حسگر، این تفاوت عملی است: کتابخانه ممکن است فقط یک نیمهٔ ماتریس را بخواند، درحالی‌که نیمهٔ دیگر از مسیرِ ناسازگار آمده باشد. این درس، مرزِ اعتبارسنجی و عامل‌گیری را با دادهٔ دوکانالهٔ تألیفی، پاسخ دستی و مرجعِ حسابِ کسری جدا می‌کند.

۱. کوواریانسِ دو کانال و معنای عامل

دو کانال را هم‌واحد و مقدارهایشان را ولتاژ فرض می‌کنیم. ماتریسِ هدف C، سطرهای چهار و دو، سپس دو و پنج دارد؛ واحدِ همهٔ خانه‌ها ولتِ مربع است. قطرها واریانس و خانه‌های خارجِ قطر کوواریانس‌اند. تعریفِ جهتِ متغیر و مشاهده و مخرجِ محاسبه در مستندات cov در NumPy 2.4 جدا آمده است؛ نامِ کوواریانس، مخرج را به‌تنهایی تعیین نمی‌کند.

عاملِ پایین‌مثلثیِ این مثال L، سطرهای دو و صفر، سپس یک و دو دارد. ضربِ L در ترانهاده‌اش C را می‌دهد: واریانسِ کانال اول چهار، کوواریانس دو، و واریانسِ دوم مجموعِ یک و چهار است. در این تفسیر، خانه‌های L واحدِ ولت دارند. اگر کانال‌ها واحدهای متفاوت داشته باشند، همین برچسبِ یکدست برای تمام خانه‌ها معتبر نیست و مدل باید دوباره تعریف شود.

در مدلِ حقیقیِ متقارن، مثبت‌معین بودن یعنی برای هر بردارِ ناصفر، فرمِ درجه‌دو مثبت باشد. توضیحِ اصلیِ نیک هایم دربارهٔ Cholesky، وجودِ تجزیه و یکتاییِ عامل با قطرِ مثبت را بیان می‌کند. عاملِ مثلثی عموماً متقارن نیست؛ آن را با ریشهٔ ماتریسیِ متقارن یا نتیجهٔ مربع‌کردنِ سادهٔ خودِ L یکی نمی‌گیریم.

ترتیبِ ضرب اهمیت دارد. در شاهدِ دستی، ضربِ ترانهادهٔ L در L، قطرهای پنج و چهار را می‌دهد و همان C نیست. آزمونِ بازسازی باید دقیقاً کمیتِ مطلوب را بسازد. مانند آزمون حساسیتِ دستگاه خطی، باید معنای معیار را پیش از مشاهدهٔ عددِ کوچک مشخص کنیم؛ اینجا دستگاهی برای مجهول‌های ولتاژ حل نمی‌کنیم.

۲. تقارنِ کامل را به حل‌گر واگذار نکنیم

قرارداد جاریِ cholesky در NumPy 2.5 صریح است: تابع، تقارنِ هرمیتی را بررسی نمی‌کند و بسته به گزینه، فقط خانه‌های یکی از مثلث‌ها و قطر را به کار می‌گیرد. برای دادهٔ حقیقیِ ما، انتخابِ پایین‌مثلثی مناسب است؛ اما انتخابِ درستِ گزینه، جای بررسیِ تقارنِ کلِ ورودی را نمی‌گیرد.

عمداً ماتریسی با سطرهای چهار و نهصدونودونه، سپس دو و پنج می‌سازیم. این ورودی متقارن نیست و خارج از پیش‌شرطِ معتبرِ تجزیه است. فراخوانیِ تشخیصیِ کتابخانه با خواندنِ نیمهٔ پایین، همان L قبلی را برمی‌گرداند؛ حاصل‌ضربِ آن C را بازسازی می‌کند، نه ماتریسِ کاملِ نامتقارن را. موفقیتِ این فراخوانی، تأییدِ تقارن یا یک API معتبر برای ورودیِ خراب نیست.

تابعِ آموزشی پیش از کتابخانه، تمام خانه‌ها را به float64 متناهی تبدیل و سپس تقارنِ عددهای تبدیل‌شده را دقیق بررسی می‌کند. هیچ میانگین‌گیریِ خاموشِ دو نیمه یا تلورانسِ پنهانی نداریم. دو عددِ اعشاری یا صحیح ممکن است هنگام تبدیل گرد شوند؛ شرطِ برابری مربوط به ماتریسِ float موجود است، نه اثباتِ برابریِ مقدارهای اندازه‌گیریِ پیش از ورود.

برای داده‌ای که بر اثر گردکردنِ مستقل، اندکی نامتقارن شده است، می‌توان سیاستِ جداگانه‌ای طراحی کرد؛ اما آستانه و معنای اصلاح باید معلوم باشد. این درس چنین سیاستی را حدس نمی‌زند. همهٔ خانه‌ها از نظر متناهی بودن نیز دیده می‌شوند، حتی اگر در مثلثِ نادیدهٔ کتابخانه باشند. ورودیِ نامعتبر، با انتخابِ نیمهٔ خوش‌رفتار معتبر نمی‌شود.

۳. کد کاملِ مرز ورودی و تجزیه

رابطِ حاضر فهرست یا تاپلِ مربعیِ غیرخالی از int یا float واقعی پایتون می‌پذیرد؛ بولی، رشته، مختلط و آرایهٔ NumPy در این رابطِ مشخص پذیرفته نیستند. این محدودیتِ ورودیِ برنامهٔ آموزشی است، نه محدودیتِ عمومیِ الگوریتم یا کتابخانه. اندازهٔ ماتریس را به دو محدود نکرده‌ایم؛ شاهدِ مهندسی دوکاناله است و مرجعِ آزمون، اندازه‌های بیشتری را نیز پوشش می‌دهد.

پس از تقارنِ کامل، عاملِ پایین‌مثلثی با قطرِ مثبت درخواست می‌شود. ماتریسِ تکین یا نامعین به مسیرِ دیگری مانند شبه‌وارون، PCA یا منظم‌سازی فرستاده نمی‌شود؛ خطای عامل‌گیری حفظ می‌شود. کوواریانسِ معتبر می‌تواند فقط مثبت‌نیمه‌معین باشد، مثلاً وقتی دو کانال کاملاً وابسته‌اند. رد شدنِ آن در این تابع، الزاماً به معنای غیرممکن بودنِ مدلِ کوواریانس نیست؛ پیش‌شرطِ این API سخت‌تر است.

"""Validate the whole real covariance before a triangle-reading Cholesky call."""
from decimal import Decimal, localcontext
from fractions import Fraction
from itertools import product
import json
import math
import platform
import numpy as np


def checked_cholesky(covariance_v2):
    """Nonempty square list/tuple of finite actual int/float; exact float symmetry."""
    if type(covariance_v2) not in (list,tuple) or not covariance_v2:
        raise ValueError('nonempty square list/tuple required')
    n = len(covariance_v2)
    converted = []
    for row in covariance_v2:
        if type(row) not in (list,tuple) or len(row) != n:
            raise ValueError('square rows required')
        values=[]
        for value in row:
            if type(value) not in (int,float):
                raise ValueError('actual int/float required, not bool')
            try:
                number=float(value)
            except OverflowError as exc:
                raise ValueError('finite representable value required') from exc
            if not math.isfinite(number):
                raise ValueError('all entries must be finite')
            values.append(number)
        converted.append(values)
    a=np.array(converted,dtype=np.float64)
    if not np.array_equal(a,a.T):
        raise ValueError('whole converted matrix must be exactly symmetric')
    # Do not average halves, add jitter, invert A, or accept semidefinite input.
    factor=np.linalg.cholesky(a,upper=False)
    if not np.isfinite(factor).all() or not np.all(np.diag(factor)>0):
        raise ArithmeticError('nonfinite or nonpositive factor diagonal')
    return factor


def exact_ldlt_oracle(values):
    """Rational unit-lower LDL^T elimination; no NumPy factorization or roots."""
    a=[[Fraction(float(v)) for v in row] for row in values]
    n=len(a); unit=[[Fraction(int(i==j)) for j in range(n)] for i in range(n)]; diagonal=[]
    for j in range(n):
        pivot=a[j][j]-sum(unit[j][k]**2*diagonal[k] for k in range(j))
        if pivot<=0: raise ValueError('not positive definite in exact rational model')
        diagonal.append(pivot)
        for i in range(j+1,n):
            unit[i][j]=(a[i][j]-sum(unit[i][k]*unit[j][k]*diagonal[k] for k in range(j)))/pivot
    return unit,diagonal


def oracle_cholesky_from_ldlt(unit,diagonal):
    with localcontext() as context:
        context.prec=100
        roots=[(Decimal(v.numerator)/Decimal(v.denominator)).sqrt() for v in diagonal]
        return np.array([[float(Decimal(v.numerator)/Decimal(v.denominator)*roots[j]) for j,v in enumerate(row)] for row in unit])


def fails(exc,values):
    try: checked_cholesky(values)
    except exc: return True
    return False


def run():
    checks={}
    def check(name,passed):
        checks[name]=bool(passed)
        assert checks[name],name
    c=((4,2),(2,5)); a=np.array(c,dtype=np.float64); factor=checked_cholesky(c)
    check('literal_factor_and_reconstruction',np.array_equal(factor,[[2,0],[1,2]]) and np.array_equal(factor@factor.T,a))
    check('multiplication_order_is_not_interchangeable',np.array_equal(factor.T@factor,[[5,2],[2,4]]) and not np.array_equal(factor.T@factor,a))
    z=np.array([[1,-1,1,-1],[1,1,-1,-1]],dtype=np.float64)
    y=factor@z
    check('literal_deterministic_coloring',np.array_equal(y,[[2,-2,2,-2],[3,1,-1,-3]]))
    check('population_zero_mean_second_moments',np.array_equal(z.mean(axis=1),[0,0]) and np.array_equal(z@z.T/4,np.eye(2)) and np.array_equal(y@y.T/4,a) and np.array_equal(np.cov(y,rowvar=True,ddof=0),a))
    check('sample_denominator_changes_target',np.allclose(np.cov(y,rowvar=True,ddof=1),a*4/3,rtol=1e-14,atol=1e-14) and not np.array_equal(np.cov(y),a))
    bad=np.array([[4,999],[2,5]],dtype=np.float64)
    diagnostic=np.linalg.cholesky(bad,upper=False)
    check('library_one_triangle_does_not_validate_full_symmetry',np.array_equal(diagnostic,factor) and not np.array_equal(diagnostic@diagnostic.T,bad) and fails(ValueError,((4,999),(2,5))))
    seen=set(); cases=0; max_factor_difference=0.; max_reconstruction_difference=0.
    for n in range(1,5):
        for off in (-1,0,1):
            b=[[int(i==j) if i<=j else off*(-1)**(i+j) for j in range(n)] for i in range(n)]
            for diag in product((1,2),repeat=n):
                matrix=tuple(tuple(sum(b[i][k]*diag[k]*b[j][k] for k in range(n)) for j in range(n)) for i in range(n))
                if matrix in seen: continue
                seen.add(matrix)
                actual=checked_cholesky(matrix)
                unit,pivots=exact_ldlt_oracle(matrix)
                reference=oracle_cholesky_from_ldlt(unit,pivots)
                # Independent exact reconstruction of the rational factorization.
                assert all(sum(unit[i][k]*pivots[k]*unit[j][k] for k in range(n))==matrix[i][j] for i in range(n) for j in range(n))
                delta=float(np.max(np.abs(actual-reference)))
                residual=float(np.max(np.abs(actual@actual.T-np.array(matrix))))
                assert np.allclose(actual,reference,rtol=2e-12,atol=2e-12)
                assert np.allclose(actual@actual.T,matrix,rtol=2e-12,atol=2e-12)
                max_factor_difference=max(max_factor_difference,delta);max_reconstruction_difference=max(max_reconstruction_difference,residual);cases+=1
    check('exact_LDLT_Decimal_root_oracle_86_unique_matrices_ONE_group',cases==86)
    check('diagonal_and_single_channel_covariances',np.array_equal(checked_cholesky(((9,),)),[[3]]) and np.array_equal(checked_cholesky(((4,0),(0,9))),[[2,0],[0,3]]))
    check('singular_semidefinite_not_accepted_as_positive_definite',fails(np.linalg.LinAlgError,((1,1),(1,1))) and fails(np.linalg.LinAlgError,((0,0),(0,0))))
    check('positive_diagonal_alone_does_not_imply_positive_definite',fails(np.linalg.LinAlgError,((1,2),(2,1))))
    check('container_shape_and_exact_numeric_type_rejections',all(fails(ValueError,v) for v in ([],[1,2],((1,2),),((True,0),(0,1)),(('1',0),(0,1)),((1+0j,0),(0,1)),np.eye(2))))
    check('all_entries_finite_even_in_unused_triangle',all(fails(ValueError,((4,v),(2,5))) for v in (float('nan'),float('inf'),-float('inf'))) and fails(ValueError,((10**1000,0),(0,1))))
    original=[[4,2],[2,5]]; result=checked_cholesky(original); result[0,0]=99
    check('input_and_next_call_not_mutated_by_returned_factor',original==[[4,2],[2,5]] and np.array_equal(checked_cholesky(original),factor))
    check('covariance_gain_four_gives_factor_gain_two',np.array_equal(checked_cholesky(((16,8),(8,20))),2*factor))
    singular=np.array([[1,1],[1,1]],dtype=np.float64); modified=singular+.25*np.eye(2); adjusted=checked_cholesky(modified.tolist())
    check('explicit_diagonal_perturbation_changes_the_matrix',np.allclose(adjusted@adjusted.T,modified,rtol=2e-12,atol=2e-12) and not np.array_equal(modified,singular))
    return {'python':platform.python_version(),'numpy':np.__version__,'units':{'covariance':'V^2','Cholesky_factor':'V','z':'dimensionless','colored_columns':'V'},'hand_covariance':a.tolist(),'hand_lower_factor':factor.tolist(),'wrong_order_LtL':(factor.T@factor).tolist(),'deterministic_Z':z.tolist(),'deterministic_Y':y.tolist(),'population_covariance_denominator_M':np.cov(y,rowvar=True,ddof=0).tolist(),'sample_covariance_denominator_M_minus_1':np.cov(y,rowvar=True,ddof=1).tolist(),'nonsymmetric_library_diagnostic_input':bad.tolist(),'library_diagnostic_factor':diagnostic.tolist(),'oracle_unique_matrices_in_ONE_group':cases,'max_factor_difference_V':max_factor_difference,'max_reconstruction_difference_V2':max_reconstruction_difference,'relative_tolerance':2e-12,'absolute_tolerance_in_corresponding_units':2e-12,'random_draws':0,'checks':checks}


if __name__=='__main__':
    print(json.dumps(run(),ensure_ascii=False,indent=2,allow_nan=False))

کد با Python 3.14.5 و NumPy 2.4.6 اجرا شده است؛ نسخهٔ مستنداتِ جاریِ عامل‌گیری 2.5 است و نسخهٔ cov خوانده‌شده 2.4. اجرای SciPy، MATLAB یا ساختِ نویزِ تصادفی نداریم. ورودی بازنویسی نمی‌شود، ولی آرایهٔ عاملِ برگشتی قابل تغییر است. تغییرِ آن، ورودی یا خروجیِ فراخوانیِ تازه را تغییر نمی‌دهد؛ ادعای تغییرناپذیریِ عمیقِ آرایه نمی‌کنیم.

۴. مرجع کسری و ساختِ شاهد دوکاناله

مرجعِ مستقل، تجزیهٔ LDLᵀ با عاملِ پایین‌مثلثیِ دارای قطرِ واحد را با Fraction انجام می‌دهد. حذفِ کسری، بدون فراخوانیِ cholesky یا محاسبهٔ ریشه در مسیرِ اصلیِ مرجع است. سپس ریشهٔ قطرهای مثبت با Decimal و دقتِ صد رقم تقریب زده می‌شود تا عاملِ قابل مقایسه ساخته شود. این مرجع نیز برای ریشهٔ گنگ تقریب دارد و گواهیِ گردشدنِ درستِ همهٔ ورودی‌ها نیست.

هشتادوشش ماتریسِ یکتای کوچک در اندازه‌های یک تا چهار ساخته و بررسی شدند؛ تمام این مجموعه یک گروهِ نام‌دار است. بازسازیِ LDLᵀ نیز با کسرهای دقیق کنترل می‌شود. بیشترین اختلافِ عامل‌ها حدود ۳٫۵۵ ضرب‌در ده به توان منفی پانزده ولت و خطای بازسازی حدود ۴٫۴۴ ضرب‌در ده به توان منفی شانزده ولتِ مربع بود. مقایسه، تلورانسِ نسبی و مطلقِ دو ضرب‌در ده به توان منفی دوازده در واحدِ مربوط دارد.

برای شاهدِ مستقلِ داده، Z دو سطر دارد: یک، منفی یک، یک، منفی یک؛ و یک، یک، منفی یک، منفی یک. ستون‌ها چهار مشاهدهٔ قطعی و بی‌بعدند. میانگینِ هر سطر صفر و حاصل‌ضربِ Z در ترانهاده‌اش تقسیم بر چهار، همانی است. این ساخت، نمونه‌گیریِ تصادفی یا اثباتِ استقلالِ آماریِ خروجیِ یک مولد نیست؛ چهار بردارِ صریح برای آزمونِ جبرند.

با Y برابر L ضرب‌در Z، سطرِ نخست دو، منفی دو، دو و منفی دو ولت و سطرِ دوم سه، یک، منفی یک و منفی سه ولت می‌شود. میانگین‌ها صفر می‌مانند و Y ضرب‌در ترانهاده‌اش با مخرج چهار دقیقاً C است. اگر cov را با پیش‌فرضِ مخرجِ تعدادمنهای‌یک صدا بزنیم، چهار‌سومِ C می‌گیریم. اختلافِ مخرج، خطای عامل‌گیری نیست.

{
  "python": "3.14.5",
  "numpy": "2.4.6",
  "units": {
    "covariance": "V^2",
    "Cholesky_factor": "V",
    "z": "dimensionless",
    "colored_columns": "V"
  },
  "hand_covariance": [
    [
      4.0,
      2.0
    ],
    [
      2.0,
      5.0
    ]
  ],
  "hand_lower_factor": [
    [
      2.0,
      0.0
    ],
    [
      1.0,
      2.0
    ]
  ],
  "wrong_order_LtL": [
    [
      5.0,
      2.0
    ],
    [
      2.0,
      4.0
    ]
  ],
  "deterministic_Z": [
    [
      1.0,
      -1.0,
      1.0,
      -1.0
    ],
    [
      1.0,
      1.0,
      -1.0,
      -1.0
    ]
  ],
  "deterministic_Y": [
    [
      2.0,
      -2.0,
      2.0,
      -2.0
    ],
    [
      3.0,
      1.0,
      -1.0,
      -3.0
    ]
  ],
  "population_covariance_denominator_M": [
    [
      4.0,
      2.0
    ],
    [
      2.0,
      5.0
    ]
  ],
  "sample_covariance_denominator_M_minus_1": [
    [
      5.333333333333333,
      2.6666666666666665
    ],
    [
      2.6666666666666665,
      6.666666666666666
    ]
  ],
  "nonsymmetric_library_diagnostic_input": [
    [
      4.0,
      999.0
    ],
    [
      2.0,
      5.0
    ]
  ],
  "library_diagnostic_factor": [
    [
      2.0,
      0.0
    ],
    [
      1.0,
      2.0
    ]
  ],
  "oracle_unique_matrices_in_ONE_group": 86,
  "max_factor_difference_V": 3.552713678800501e-15,
  "max_reconstruction_difference_V2": 4.440892098500626e-16,
  "relative_tolerance": 2e-12,
  "absolute_tolerance_in_corresponding_units": 2e-12,
  "random_draws": 0,
  "checks": {
    "literal_factor_and_reconstruction": true,
    "multiplication_order_is_not_interchangeable": true,
    "literal_deterministic_coloring": true,
    "population_zero_mean_second_moments": true,
    "sample_denominator_changes_target": true,
    "library_one_triangle_does_not_validate_full_symmetry": true,
    "exact_LDLT_Decimal_root_oracle_86_unique_matrices_ONE_group": true,
    "diagonal_and_single_channel_covariances": true,
    "singular_semidefinite_not_accepted_as_positive_definite": true,
    "positive_diagonal_alone_does_not_imply_positive_definite": true,
    "container_shape_and_exact_numeric_type_rejections": true,
    "all_entries_finite_even_in_unused_triangle": true,
    "input_and_next_call_not_mutated_by_returned_factor": true,
    "covariance_gain_four_gives_factor_gain_two": true,
    "explicit_diagonal_perturbation_changes_the_matrix": true
  }
}

پانزده گروه موفق شدند: عاملِ دستی، ترتیب ضرب، شاهدِ ستون‌ها، مخرج‌ها، دامِ یک نیمه، مرجعِ کسری و قراردادهای خطا پوشش دارند. آزمون پروژکتور PCA جهتِ بیشترین واریانس را انتخاب می‌کند؛ اینجا تمام کوواریانسِ هدف عامل‌گیری می‌شود و حذفِ مؤلفه یا بهبودِ دسته‌بندی نداریم. تلورانسِ کوچک، عدم قطعیتِ دادهٔ واقعی یا مقاومتِ همهٔ ماتریس‌های نزدیک به تکینی را تعیین نمی‌کند.

۵. چهار پرسش و مرز مثبت‌معین بودن

پرسش اول: اگر C چهار برابر شود، عامل چه می‌شود؟ پاسخ: دو برابرِ L با قطرِ مثبت، عاملِ همان ماتریسِ تازه است، چون حاصل‌ضرب دو عامل چهار برابر می‌شود. این تغییر، با دو برابر شدنِ دامنهٔ هر دو کانال سازگار است. عوض‌کردنِ تنها یک کانال، همهٔ خانه‌های کوواریانس را یکسان مقیاس نمی‌کند.

پرسش دوم: قطرِ مثبت برای موفقیت کافی است؟ پاسخ: خیر؛ ماتریسِ متقارن با قطرهای یک و خانه‌های خارجِ قطرِ دو، در جهتِ یک و منفی یک فرمِ درجه‌دوی منفی می‌دهد. در این شاهد مقدار منفی دو است و عامل‌گیری رد می‌شود. شکلِ مربع و واریانس‌های مثبت، به‌تنهایی وابستگیِ سازگارِ دو کانال را ثابت نمی‌کنند.

پرسش سوم: چرا cov پیش‌فرض با C برابر نشد؟ پاسخ: چهار مشاهدهٔ شاهد را جمعیتِ محدودِ همین تمرین با مخرج چهار گرفته‌ایم؛ پیش‌فرض، مخرج سه دارد. در نتیجه همهٔ خانه‌ها چهار‌سوم می‌شوند. انتخابِ مخرج باید تابعِ هدفِ تحلیل باشد؛ تغییرِ آن برای شبیه شدنِ تصادفیِ خروجی‌ها، معیارِ آزمونِ روشنی نمی‌سازد.

پرسش چهارم: می‌توان برای رفعِ خطا، خودکار عددی به قطر افزود؟ پاسخ: این تابع چنین نمی‌کند. افزودنِ یک‌چهارمِ همانی به شاهدِ تکین، ماتریسِ دیگری با واریانس‌های متفاوت می‌سازد؛ برنامه همین تغییرِ صریح را جدا نشان می‌دهد. این کار گواهیِ نویزِ فیزیکی یا اصلاحِ بی‌اثر نیست. منابع در ۲۶ مهر ۱۴۰۵ خوانده شده‌اند؛ نتیجه، آزمونِ عددیِ تألیفی است و تصاویر مفهومی‌اند.

منابع و مطالعهٔ بیشتر

• NumPy؛ قرارداد Cholesky و خواندن یک نیمه
• نیک هایم؛ تعریف و یکتایی تجزیهٔ Cholesky
• NumPy؛ جهتِ مشاهده‌ها و مخرجِ کوواریانس

بررسی منابع: ۲۶ مهر ۱۴۰۵. تحلیل‌ها و پیشنهادهای اجرایی، جمع‌بندی نویسنده‌اند.