QR分解失效时的求根问题解决方法咨询

Great question! The issue you're hitting with companion matrix-based QR decomposition for polynomial root-finding is a common one tied to the numerical behavior of both companion matrices and basic QR iteration.

Companion matrices are highly non-normal, meaning their condition numbers are often very large—this makes basic QR iteration (which converges slowly for non-normal matrices) stall or fail to converge properly, even for polynomials with simple real roots like your examples. For polynomials with paired roots (like ±1, ±3 in your first equation), basic QR struggles to separate these closely grouped eigenvalues, leading to the "unsolvable" behavior you’re seeing.

If you want to stick with the QR-based method, here are targeted modifications to make it work for your cases:

1. Balance the Companion Matrix First #

Before running QR iteration, apply a balancing transformation to reduce the companion matrix’s condition number. The goal is to find a diagonal scaling matrix D such that D⁻¹ * C * D (where C is the companion matrix) has balanced row and column L1 norms. This simple preprocessing step drastically improves numerical stability and convergence speed, especially for polynomials with roots of varying magnitudes.

2. Use Shifted QR Iteration with Wilkinson Shifts #

Basic QR iteration converges linearly, which is too slow for many cases. Switching to Wilkinson-shifted QR iteration boosts convergence to quadratic speed for simple real roots, and handles complex conjugate root pairs gracefully. For your examples, this shifted approach will quickly converge to all four real roots without stalling. When the iteration detects a 2x2 subblock corresponding to complex conjugate roots, you can use implicit complex shifts or block QR iteration to handle these subblocks directly (avoiding unnecessary complex arithmetic in real-valued implementations).

3. Factor the Polynomial to Remove Clustered/Repeated Roots #

While your examples don’t have repeated roots, this step is critical for cases that do. Compute the greatest common divisor (GCD) of the polynomial and its derivative to split the original polynomial into factor polynomials with no repeated roots. Then run QR iteration on the companion matrix of each factor separately. For example, if you had a polynomial like (x-2)²(x+3), the GCD would isolate x-2, letting you split the problem into solving x-2 and (x-2)(x+3) independently.

If modifying the QR approach feels cumbersome, these specialized polynomial root-finding algorithms are more reliable for all cases, including yours:

Jenkins-Traub Algorithm #

This is the industry standard for numerical polynomial root-finding, optimized for stability and efficiency. It uses three phases:

• Initialization: Builds a linear iteration sequence to prepare for root separation
• Shift phase: Uses fixed shifts to isolate the smallest-magnitude root
• Convergence phase: Switches to variable shifts for fast, stable convergence

It handles real roots, complex roots, repeated roots, and clustered roots seamlessly—perfect for your two example polynomials.

Laguerre's Method #

A single-root iterative method with cubic convergence (faster than Newton’s method). After finding one root, you reduce the polynomial’s degree via polynomial division and iterate for the next root. It works well for both real and complex roots, and is less sensitive to initial guesses than Newton’s method.

Durand-Kerner Method #

A simultaneous iteration method that estimates all roots at once, no polynomial division required. You initialize guesses for all roots, then iteratively update each guess until convergence. It’s easy to implement and performs well for polynomials with multiple closely grouped roots.

For your first polynomial x⁴ -10x² +9, any of these modified or alternative methods will quickly return the roots 1, -1, 3, -3. For the second polynomial x³ -5x² -4x +20, they’ll converge to 2, -2, 5 without issues.

内容的提问来源于stack exchange，提问作者FROST SKIPPER