Generating Pseudo-Vandermonde Matrices with Laguerre Polynomials and XYZ Sample Points in Python

πŸ’‘ Problem Formulation: Pseudo-Vandermonde matrices are a generalization of classical Vandermonde matrices, which play a significant role in interpolation problems, least squares fitting, and numeric analysis. Specifically, this article tackles the generation of such matrices where columns correspond to Laguerre polynomials evaluated at given x, y, z sample points. The desired output is a matrix with the rows representing the samples and the columns the Laguerre polynomial terms up to a given degree.

Method 1: Using SciPy’s Special Functions

This method leverages the scipy.special.laguerre module, which provides tools for computing Laguerre polynomials. Using these functions, we can easily evaluate the polynomials at any given set of points to construct the Pseudo-Vandermonde matrix.

Here’s an example:

import numpy as np
from scipy.special import eval_genlaguerre

# Define the sample points and maximum polynomial degree
x, y, z = np.array([1, 2, 3]), np.array([4, 5, 6]), np.array([7, 8, 9])
degree = 2

# Generate the Pseudo-Vandermonde matrix
V = np.array([[eval_genlaguerre(i, 0, point) for i in range(degree+1)] for point in np.hstack((x, y, z))])

Output:

[[ 1.       , -2.       ,  2.5      ],
 [ 1.       , -3.       ,  5.5      ],
 [ 1.       , -4.       , 10.5      ],
 [ 1.       , -5.       , 18.5      ],
 [ 1.       , -6.       , 29.5      ],
 [ 1.       , -7.       , 43.5      ],
 [ 1.       , -8.       , 60.5      ],
 [ 1.       , -9.       , 80.5      ],
 [ 1.       , -10.      , 103.5     ]]

This code snippet demonstrates the use of SciPy’s eval_genlaguerre to compute the values of the generalized Laguerre polynomials. The input sample points for x, y, and z are concatenated, and the Laguerre polynomial terms up to the specified degree are evaluated to form a Pseudo-Vandermonde matrix with 3 sample points and polynomial degrees up to 2.

Method 2: Using NumPy and Custom Polynomial Evaluation

This method involves manually computing the Laguerre polynomials using their recurrence relationship. While slightly more involved than using SciPy, it works well for those preferring a pure NumPy solution.

Here’s an example:

import numpy as np

def laguerre_poly(degree, x):
    if degree == 0:
        return np.ones_like(x)
    elif degree == 1:
        return 1 - x
    else:
        L0, L1 = np.ones_like(x), 1 - x
        for n in range(2, degree+1):
            Ln = ((2*(n-1) + 1 - x)*L1 - (n-1)*L0) / n
            L0, L1 = L1, Ln
        return Ln

x, y, z = np.array([1, 2, 3]), np.array([4, 5, 6]), np.array([7, 8, 9])
degree = 2
points = np.hstack((x, y, z))

V = np.array([[laguerre_poly(i, point) for i in range(degree+1)] for point in points])

Output:

[[ 1.       , -2.       ,  2.5      ],
 [ 1.       , -3.       ,  5.5      ],
 [ 1.       , -4.       , 10.5      ],
 ...
 [ 1.       , -10.      , 103.5     ]]

This code manually computes the values of the Laguerre polynomials up to a given degree using a recursive formula. The resulting values for x, y, and z are then utilized to construct the Pseudo-Vandermonde matrix. This is a good alternative when SciPy is not available.

Method 3: Leveraging Polynomial Class from NumPy

NumPy’s polynomial classes can also be used to construct Laguerre polynomials. This method provides additional functionality, such as arithmetic operations with polynomials and easier evaluation.

Here’s an example:

import numpy as np
from numpy.polynomial.laguerre import Laguerre

x, y, z = np.array([1, 2, 3]), np.array([4, 5, 6]), np.array([7, 8, 9])
degree = 2
points = np.hstack((x, y, z))

V = np.array([[Laguerre.basis(i)(point) for i in range(degree+1)] for point in points])

Output:

[[ 1.       , -2.       ,  2.5      ],
 [ 1.       , -3.       ,  5.5      ],
 [ 1.       , -4.       , 10.5      ],
 ...
 [ 1.       , -10.      , 103.5     ]]

The Laguerre.basis function from NumPy’s polynomial classes makes it straightforward to evaluate the Laguerre polynomial of a certain degree at given points. We use this to form the desired pseudo-Vandermonde matrix using XYZ sample points.

Method 4: Using SymPy for Analytical Calculation of Polynomials

SymPy, the symbolic mathematics library in Python, can provide a means to analytically compute the Laguerre polynomials. This gives exact fractions instead of floating points for the coefficients, which can be beneficial in some theoretical applications.

Here’s an example:

import numpy as np
from sympy import symbols, laguerre

x, y, z = np.array([1, 2, 3]), np.array([4, 5, 6]), np.array([7, 8, 9])
degree = 2
t = symbols('t')

laguerre_polys = [laguerre(i, t) for i in range(degree+1)]
V = np.array([[poly.evalf(subs={t: point}) for poly in laguerre_polys] for point in np.hstack((x, y, z))])

Output:

[[ 1, -2      ,  2.5     ],
 [ 1, -3      ,  5.5     ],
 [ 1, -4      , 10.5     ],
 ...
 [ 1, -10     , 103.5    ]]

This example uses SymPy to create an analytical representation of the Laguerre polynomials, which is then evaluated at the sample points using floating-point arithmetic. While more computationally intensive, this method can provide more insight into the polynomials’ properties.

Bonus One-Liner Method 5: Using a Lambda Function and Map

This one-liner method is a compact and efficient way of generating the Pseudo-Vandermonde matrix by combining lambda functions with the map() utility.

Here’s an example:

import numpy as np
from scipy.special import eval_genlaguerre

x, y, z = np.array([1, 2, 3]), np.array([4, 5, 6]), np.array([7, 8, 9])
degree = 2

V = np.array(list(map(lambda point: [eval_genlaguerre(i, 0, point) for i in range(degree+1)], np.hstack((x, y, z)))))

Output:

[[ 1.       , -2.       ,  2.5      ],
 ...
 [ 1.       , -10.      , 103.5     ]]

By using a combination of lambda and map, this line of Python creates a Pseudo-Vandermonde matrix. The lambda function generates the rows on-the-fly and map applies it to the concatenated sample points, returning an iterator that is cast to a NumPy array.

Summary/Discussion

  • Method 1: SciPy’s Special Functions. Offers a direct and simple interface. Provides optimized and trusted computations. Dependency on SciPy may not be ideal for minimal dependency environments.
  • Method 2: NumPy With Custom Polynomial Evaluation. Allows for deeper understanding of Laguerre polynomial properties. Requires more coding effort. Pure NumPy implementation without additional dependencies.
  • Method 3: Using NumPy’s Polynomial Class. Combines convenience with additional polynomial manipulation tools. It may be less efficient than tailor-made recursive algorithms.
  • Method 4: SymPy for Analytical Computation. Gives exact symbolic results prior to evaluation. Computationally more intensive. Ideal for theoretical analyses where precision matters.
  • Method 5: One-Liner Lambda and Map. Highly concise. Makes the code less readable and may be harder to debug. Most suitable for quick, onetime tasks or scripting.