ब्लॉक GMRES क्यों अटकता है: रेज़िडुअल ब्लॉक का QR ही सही उपाय

ब्लॉक GMRES में कई RHS पर ठहराव/डाइवर्जेंस का सही कारण और QR‑आधारित रेज़िडुअल से इसे कैसे ठीक करें। सही लीस्ट‑स्क्वेयर्स रिस्टार्ट, संक्षिप्त कोड और तुलना।

ब्लॉक GMRES दिखने में सिंगल‑RHS GMRES से बहुत मिलता‑जुलता लगता है: वही Krylov उपस्थान का विचार, वही लीस्ट‑स्क्वेयर्स चरण, और मिलती‑जुलती रिस्टार्ट रणनीति। पेच इस बात में है कि शुरुआती रेज़िडुअल को छोटे लीस्ट‑स्क्वेयर्स प्रॉब्लम में कैसे समाहित किया जाता है। यह हिस्सा अगर ठीक न हो, तो कई दाएँ‑हाथी पक्षों के लिए विधि ठहर सकती है या यहाँ तक कि भटक भी सकती है। नीचे एक न्यूनतम, काम करने वाला वॉकथ्रू है, जो गलती और उसका सटीक सुधार दिखाता है।

समस्या की सेटिंग और एक न्यूनतम GMRES बेसलाइन

हम संदर्भ के लिए एक सिंगल‑RHS, कॉम्पैक्ट रिस्टार्टेड GMRES से शुरुआत करेंगे। यह मानक विवरणों के अनुसार Arnoldi प्रक्रिया और एक छोटे लीस्ट‑स्क्वेयर्स हल का उपयोग करता है।

import numpy as np
def gmres_cycle(M, rhs, x_init=None, eps=1e-5, max_steps=100, krylov_dim=30):
    n = len(rhs)
    sol = np.zeros(n) if x_init is None else x_init.copy()
    res = rhs - M @ sol
    beta0 = np.linalg.norm(res)
    if beta0 < eps:
        return sol
    for outer_k in range(0, max_steps, krylov_dim):
        basis = np.zeros((n, krylov_dim + 1))
        hess = np.zeros((krylov_dim + 1, krylov_dim))
        basis[:, 0] = res / beta0
        for j in range(krylov_dim):
            wv = M @ basis[:, j]
            for i in range(j + 1):
                hess[i, j] = np.dot(basis[:, i], wv)
                wv -= hess[i, j] * basis[:, i]
            hess[j + 1, j] = np.linalg.norm(wv)
            if hess[j + 1, j] < 1e-14:
                break
            basis[:, j + 1] = wv / hess[j + 1, j]
        g = np.zeros(j + 2)
        g[0] = beta0
        Hsub = hess[:j + 2, :j + 1]
        coef, *_ = np.linalg.lstsq(Hsub, g, rcond=None)
        sol += basis[:, :j + 1] @ coef
        res = rhs - M @ sol
        beta0 = np.linalg.norm(res)
        if beta0 < eps:
            return sol
    return sol

साधारण ब्लॉक GMRES और गलती कहाँ होती है

सीधे “पोर्ट” में, वेक्टर रेज़िडुअल की जगह रेज़िडुअल ब्लॉक लिया जाता है और Arnoldi की जगह ब्लॉक Arnoldi। लीस्ट‑स्क्वेयर्स को Frobenius नॉर्म में हल किया जाता है। नीचे दिए गए दो हिस्से एक ब्लॉक Arnoldi रूटीन और एक रिस्टार्टेड ब्लॉक GMRES परिभाषित करते हैं, जो सिंगल‑RHS संरचना का प्रतिबिंब है।

def arnoldi_block(M, W0, mstep):
    n, blk = W0.shape
    W = np.zeros((n, blk * (mstep + 1)))
    T = np.zeros((blk * (mstep + 1), blk * mstep))
    Q, _ = np.linalg.qr(W0, mode='reduced')
    W[:, :blk] = Q
    for j in range(mstep):
        Bj = W[:, j * blk : (j + 1) * blk]
        Zj = M @ Bj
        for i in range(j + 1):
            Bi = W[:, i * blk : (i + 1) * blk]
            Tij = Bi.T @ Zj
            T[i * blk : (i + 1) * blk, j * blk : (j + 1) * blk] = Tij
            Zj -= Bi @ Tij
        Qj, Rj = np.linalg.qr(Zj, mode='reduced')
        W[:, (j + 1) * blk : (j + 2) * blk] = Qj
        T[(j + 1) * blk : (j + 2) * blk, j * blk : (j + 1) * blk] = Rj
    return W, T
def gmres_block_cycle(M, RHS, X_init=None, eps=1e-5, max_steps=100, krylov_dim=30):
    n, blk = RHS.shape
    X = np.zeros((n, blk)) if X_init is None else X_init.copy()
    RES = RHS - M @ X
    if np.linalg.norm(RES, ord='fro') < eps:
        return X
    for outer_k in range(0, max_steps, krylov_dim):
        W, T = arnoldi_block(M, RES, krylov_dim)
        G = np.zeros(((krylov_dim + 1) * blk, blk))
        G[:blk, :] = np.linalg.norm(RES, axis=0).reshape(-1, 1) * np.eye(blk)
        Tsub = T[:(krylov_dim + 1) * blk, :krylov_dim * blk]
        Ycoef, *_ = np.linalg.lstsq(Tsub, G, rcond=None)
        X += W[:, :krylov_dim * blk] @ Ycoef
        RES = RHS - M @ X
        if np.linalg.norm(RES, ord='fro') < eps:
            return X
    return X

मुद्दा सूक्ष्म है और अपडेट साइन पलटने से इसका कोई संबंध नहीं। सिंगल‑RHS मामले में, लीस्ट‑स्क्वेयर्स के दाएँ‑हाथी पक्ष में पहला कैनोनिकल बेसिस वेक्टर स्केलर β से गुणा होता है। इसे गलत तरीके से सामान्यीकृत कर देना—यानी रेज़िडुअल ब्लॉक के कॉलम‑वाइस नॉर्म से बनी एक डायगोनल रख देना—छोटे सिस्टम को असंगत बना देता है। तब ऑप्टिमाइज़र गलत मात्रा को न्यूनतम करता है, और यही समझाता है कि एक से अधिक दाएँ‑हाथी पक्ष होते ही रेज़िडुअल बढ़ सकता है या ठहर सकता है।

मुख्य बात: रेज़िडुअल ब्लॉक का QR ही छोटे प्रॉब्लम को संचालित करता है

ब्लॉक GMRES में वर्तमान रेज़िडुअल ब्लॉक का QR फैक्टराइज़ेशन जरूरी है। Q0, R0 = qr(R) निकालें, Arnoldi प्रक्रिया के लिए Q0 को शुरुआती ब्लॉक के रूप में लें, और लीस्ट‑स्क्वेयर्स के दाएँ‑हाथी पक्ष के ऊपर‑बाएँ ब्लॉक में R0 को सीधे रखें। यही सिंगल‑RHS मामले से मूलभूत अंतर है, जहाँ स्केलर β वही भूमिका निभाता है।

सही लीस्ट‑स्क्वेयर्स सेटअप के साथ सुधरा हुआ ब्लॉक GMRES

नीचे दिया गया संस्करण ठीक यही करता है। बाकी एल्गोरिथ्म अपरिवर्तित रहता है।

def gmres_block_cycle_fixed(M, RHS, X_init=None, eps=1e-5, max_steps=100, krylov_dim=30):
    n, blk = RHS.shape
    X = np.zeros((n, blk)) if X_init is None else X_init.copy()
    RES = RHS - M @ X
    if np.linalg.norm(RES, ord='fro') < eps:
        return X
    for outer_k in range(0, max_steps, krylov_dim):
        Q0, R0 = np.linalg.qr(RES, mode='reduced')
        W, T = arnoldi_block(M, Q0, krylov_dim)
        G = np.zeros(((krylov_dim + 1) * blk, blk))
        G[:blk, :blk] = R0
        Tsub = T[:(krylov_dim + 1) * blk, :krylov_dim * blk]
        Ycoef, *_ = np.linalg.lstsq(Tsub, G, rcond=None)
        X += W[:, :krylov_dim * blk] @ Ycoef
        RES = RHS - M @ X
        if np.linalg.norm(RES, ord='fro') < eps:
            return X
    return X

पुनरुत्पाद्य परीक्षण ढांचा

यह छोटा सा सेटअप ऊपर चर्चा की गई स्थिति को प्रतिबिंबित करता है। यह त्वरित जाँच के लिए और उसी मैट्रिक्स पर सिंगल‑RHS तथा ब्लॉक रूपों की तुलना करने में उपयोगी है।

np.random.seed(42)
n = 20
A = np.diag(np.linspace(1, 100, n)) + 0.001 * np.random.randn(n, n)
A = (A + A.T) / 2
b = np.random.randn(n, 1)
x_ref = gmres_cycle(A, b[:, 0], krylov_dim=10)
X_blk = gmres_block_cycle_fixed(A, b, krylov_dim=10)

यह बारीकी क्यों मायने रखती है

ब्लॉक रूप Frobenius नॉर्म को ब्लॉक Krylov उपस्थान पर न्यूनतम करता है। रिस्टार्ट के समय छोटे लीस्ट‑स्क्वेयर्स सिस्टम में सही रेज़िडुअल ज्यामिति का एन्कोड होना चाहिए। रेज़िडुअल ब्लॉक का QR, R0 के माध्यम से, वही ज्यामिति कैप्चर करता है, और छोटे सिस्टम के दाएँ‑हाथी पक्ष में R0 पहुँचाने से मिनिमाइज़ेशन वास्तविक उद्देश्य के साथ संरेखित रहता है। डायगोनल नॉर्म रखना देखने में निर्दोष लगता है, लेकिन यह संरेखण तोड़ देता है।

निष्कर्ष

अगर ब्लॉक GMRES कई दाएँ‑हाथी पक्षों पर असफल हो रहा है जबकि सिंगल‑RHS संस्करण ठीक काम करता है, तो रिस्टार्ट पर छोटे लीस्ट‑स्क्वेयर्स सिस्टम के निर्माण की जाँच करें। Q0, R0 = qr(R) निकालें, Arnoldi को Q0 से सीड करें, और दाएँ‑हाथी पक्ष के अग्रणी ब्लॉक में R0 रखें। ऐसा करते ही अपडेट की दिशा पहले से सही होती है और साइन पलटने की कोई जरूरत नहीं। Python और C++ इम्प्लीमेंटेशन की तुलना करते समय, मैट्रिक्स ऑपरेशन्स को ऊपर दिए कोड की तरह ही रखें; यह विधि डिजाइन से iterative है, और Arnoldi चरणों के लिए स्पष्ट लूप पूरी तरह उपयुक्त हैं।