Conjugate Residual Method in Polynomial Space

Classes
class	PolynomialFilterInterface
	An interface to define R(S), i.e. set R and S, and compute R(S)*v, where R is the residual polynomial, S is a sparse matrix, and v is a vector. More...

Functions
Matrix	FilteredConjugateResidualPolynomial (const Matrix &pp, const Vector &intv, const Vector &w, mkIndex niter)
	This routine employs a conjugate-residual-type algorithm in polynomial space to minimize ||P(z)-Q(z)||w, where P(z), the base filter, is the input piecewise polynomial, and Q(z) is the output polynomial satisfying Q(0)==1, i.e. the constant term of Q(z) is 1. More...
Vector	FilteredConjugateResidualMatrixPolynomialVectorProduct (const SparseMatrix &A, const Vector &x0, const Vector &b, const Matrix &pp, const Vector &intv, const Vector &w, mkIndex niter, Real tol=0.0)
	This routine employs a conjugate-residual-type algorithm in polynomial space to compute x=x0+s(A)*r0 with r0=b-A*x0, such that ||1-z*s(z)-P(z)||w is minimized. More...

Detailed Description

Function Documentation

FilteredConjugateResidualMatrixPolynomialVectorProduct()

Vector FilteredConjugateResidualMatrixPolynomialVectorProduct	(	const SparseMatrix &	A,
		const Vector &	x0,
		const Vector &	b,
		const Matrix &	pp,
		const Vector &	intv,
		const Vector &	w,
		mkIndex	niter,
		Real	tol = 0.0
	)

This routine employs a conjugate-residual-type algorithm in polynomial space to compute x=x0+s(A)*r0 with r0=b-A*x0, such that ||1-z*s(z)-P(z)||_w is minimized.

• P(z) is a given piecewise polynomial, called the base filter.
  P(z) is expanded in the ‘translated’ (scale-and-shift) Chebyshev basis for each interval, and presented as a matrix of Chebyshev coefficients pp.
• s(z) is a polynomial of degree up to niter, the number of conjugate-residual iterations.

Parameters

A	is a sparse matrix.
x0,b	are vectors.
niter	is the number of conjugate-residual iterations. Therefore, the degree of Q(z) is up to niter+1.
intv	is a vector which defines the intervals. The jth interval is [intv(j),intv(j+1)).
w	is a vector of Chebyshev weights. The weight of jth interval is w(j). The interval weights define the inner product of two continuous functions and then the derived w-norm ||P(z)-Q(z)||w.
pp	is a matrix of Chebyshev coefficients which defines the piecewise polynomial P(z). For z in the jth interval [intv(j), intv(j+1)), P(z) equals Pj(z) = pp(1,j)*S0(z) + pp(2,j)*S1(z) + ... + pp(niter+2,j)*Sniter+1(z),
tol	is the tolerance; if the residual polynomial in z-norm is dropped by a factor lower than tol, then stop the conjugate-residual iteration.

Remarks

Here S_i(z) is the ‘translated’ Chebyshev polynomial in the jth interval, with S_i((z-c)/h) = T_i(z), c = (intv(j) + intv(j+1))/2, h = (intv(j+1) - intv(j))/2,
where T_i(z) is the Chebyshev polynomial of the first kind, T₀(z) = 1, T₁(z) = z, and T_i(z) = 2*z*T_i-1(z) - T_i-2(Z) for i>=2.

Returns

A vector x=x0+s(A)*r0 with r0=b-A*x0, such that ||(1-P(z))-z*s(z)||_w is minimized, subject to that s(z) is a polynomial of degree up to niter, where P(z) is the base filter. In short, z*s(z) approximates 1-P(z).

Remarks

Since z*s(z) approximates 1-P(z), P(0)==1 is expected. If P(0) 1, one can translate P(z). For example, if P(0)==0, one can use 1-P(z) as input instead of P(z).

Typically, the base filter, defined by pp and intv, is from Hermite interpolation in intervals [intv(j),intv(j+1)) for j=1,...,nintv, with nintv the number of intervals.

A common application is to compute R(A)*b, where R(z) approximates 1-P(z). In this case, one can set x0 = 0 and then the return vector is x = s(A)*b, where z*s(z) approximates 1-P(z); therefore, A*x is the wanted R(A)*b.

See also

PolynomialFilterInterface::filteredSparseMatrixPolynomialVectorProduct().

FilteredConjugateResidualPolynomial()

Matrix FilteredConjugateResidualPolynomial	(	const Matrix &	pp,
		const Vector &	intv,
		const Vector &	w,
		mkIndex	niter
	)

This routine employs a conjugate-residual-type algorithm in polynomial space to minimize ||P(z)-Q(z)||_w, where P(z), the base filter, is the input piecewise polynomial, and Q(z) is the output polynomial satisfying Q(0)==1, i.e. the constant term of Q(z) is 1.

Both P(z) and Q(z) are expanded in the ‘translated’ (scale-and-shift) Chebyshev basis for each interval and presented as matrices of Chebyshev coefficients, denoted by pp and qq, respectively.

Parameters

niter	is the number of conjugate-residual iterations. Therefore, the degree of Q(z) is up to niter+1.
intv	is a vector which defines the intervals. The jth interval is [intv(j),intv(j+1)).
w	is a vector of Chebyshev weights. The weight of jth interval is w(j). The interval weights define the inner product of two continuous functions and then the derived w-norm ||P(z)-Q(z)||w.
pp	is a matrix of Chebyshev coefficients which defines the piecewise polynomial P(z). For z in the jth interval [intv(j), intv(j+1)), P(z) equals Pj(z) = pp(1,j)*S0(z) + pp(2,j)*S1(z) + ... + pp(niter+2,j)*Sniter+1(z),

Remarks

Here S_i(z) is the ‘translated’ Chebyshev polynomial in the jth interval, with S_i((z-c)/h) = T_i(z), c = (intv(j) + intv(j+1))/2, h = (intv(j+1) - intv(j))/2,
where T_i(z) is the Chebyshev polynomial of the first kind, T₀(z) = 1, T₁(z) = z, and T_i(z) = 2*z*T_i-1(z) - T_i-2(Z) for i>=2.

Returns

A matrix, denoted by qq, represents a polynomial Q(z) with degree up to 1+niter and satisfying Q(0)==1, such that ||P(z))-Q(z)||_w is minimized.
This polynomial Q(z) is expanded in the ‘translated’ Chebyshev basis for each interval.
To be precise, considering z in [intv(j),intv(j+1)), Q(z) equals Q_j(z) = qq(1,j)*S₀(z) + qq(2,j)*S₁(z) + ... + qq(niter+2,j)*S_niter+1(z).

Remarks

Since Q(0)==1, P(0)==1 is expected. If P(0) 1, one can translate P(z). For example, if P(0)==0, one can use 1-P(z) as input instead of P(z).

Typically, the base filter, defined by pp and intv, is from Hermite interpolation in intervals [intv(j),intv(j+1)) for j=1,...,nintv, with nintv the number of intervals.