3
$\begingroup$

Recently I have posted a question on SO, but maybe here is better place to ask. So, I have data and I want to fit polynomial of order $k$ orthogonal with respect to weighted inner product of functions: $$ \left<f, g\right> = \sum_{j}^{n} \omega_j f(x_j) g(x_j) $$ Here $w_j$ is a vector of weights (constants assigned to each of $n$ points $x_j$).

When all $w_j = 1$ we have ordinary euclidean inner product $\left<f, g\right> = \sum f(x_j) g(x_j)$ and we can do it by QR decomposition of Vandermonde matrix $V_{ij}$:

$$ V = \begin{bmatrix} 1 & x_1 & x_1^2 &... & x_1^k \\ 1 & x_2 & x_2^2 &... & x_2^k\\ ... & ... & ... & ...\\ 1 & x_n & x_n^2 & ... & x_n^k \end{bmatrix}.$$ QR decomposition: $$ V = QR, \quad QQ^T = 1. $$ So the columns of $Q$ are orthogonal to each other with respect to ordinary euclidean inner product $\left<f, g\right> = \sum f(x_j) g(x_j)$.

My question is: How can I use QR decomposition to get the matrix, whose columns are orthogonal with respect to weighted inner product of functions?

I was thinking about introducing diagonal weights matrix

$$ W = \begin{bmatrix} \omega_1 & 0 & ...\\ 0 & \omega_2 & ...\\ ... & ... & ... \end{bmatrix}$$ and applying QR decomposition to matrix $\sqrt{W}V$:

$$ \sqrt{W}V = \widetilde{Q}R, \quad \widetilde{Q}\widetilde{Q}^T = 1 $$ and after that defining Q as $$ Q = \left(\sqrt{W}\right)^{-1}\widetilde{Q}, $$ but I haven't succeeded proving that and maybe it's wrong.

Could anyone help me? I know how to generate such polynomial recursively, but I want to know if one could somehow use QR decomposition.

$\endgroup$
3
$\begingroup$

Not sure about the QR decomposition (see edit), but I found a way to do it with an eigendecomposition (derivation below).

Formulating the problem

Let $n \times (k+1)$ matrix $V$ be the Vandermonde matrix, as defined in the original question. We want to factor $V$ into a product of two matrices: 1) An $n \times (k+1)$ matrix $A$, where $A_{ij}$ contains the value of the $j$th orthogonal polynomial evaluated at the $i$th data point, and 2) A $(k+1) \times (k+1)$ matrix $B$:

$$A B = V$$

The columns of $A$ must be orthogonal in the sense defined in the original question. That is:

$$A^T W A = I$$

where $I$ is the identity matrix and $W$ is a diagonal $n \times n$ matrix containing the weights. The solution below actually applies for any symmetric, positive definite $W$; it doesn't have to be diagonal.

Solution

One possible solution is:

$$A = V U \Lambda^{-\frac{1}{2}} \quad B = \Lambda^{\frac{1}{2}} U^T$$

where $U \Lambda U^T$ is the eigendecomposition of $V^T W V$

The solution is actually only specified up to unitary transformations (i.e. rotations and sign flips), so there are infinitely many. That is, given a solution $(A,B)$ and any matrix $M$ where $M^T M = I$, then $(A M, M^T B)$ is also a solution.

Derivation

Since $AB = V$ we can write:

$$(A B)^T W (A B) \enspace = \enspace B^T A^T W A B \enspace = \enspace V^T W V$$

Since $A^T W A = I$ this simplifies to:

$$B^T B = V^T W V$$

Let $U$ be an orthogonal matrix containing the eigenvectors of $V^T W V$ on the columns, and let diagonal matrix $\Lambda$ contain the corresponding eigenvalues. So $U \Lambda U^T = V^T W V$. Then a solution to the above equation is:

$$B = \Lambda^\frac{1}{2} U^T$$

Plugging this back into $A B = V$ gives:

$$A \Lambda^{\frac{1}{2}} U^T = V$$

Since $U$ is orthogonal and $\Lambda$ is diagonal, we can solve for $A$ by right-multiplying both sides by $U \Lambda^{-\frac{1}{2}}$:

$$A = V U \Lambda^{-\frac{1}{2}}$$

Edit: Solution based on QR decomposition

The OP proposed the following solution (I've changed some variable names): Let $\tilde{Q} \tilde{R}$ be the QR decomposition of $W^{\frac{1}{2}} V$. Then:

$$A = W^{-\frac{1}{2}} \tilde{Q}$$

Here's a proof that the columns are orthogonal:

$$A^T W A \tag{This must equal the identity matrix}$$

$$= \tilde{Q}^T (W^{-\frac{1}{2}})^T W W^{-\frac{1}{2}} \tilde{Q} \tag{Substitute in expression for $A$}$$

$$= \tilde{Q}^T W^{-\frac{1}{2}} W W^{-\frac{1}{2}} \tilde{Q} \tag{$W$ is symmetric}$$

$$= \tilde{Q}^T \tilde{Q} \tag{$W^{-\frac{1}{2}} W W^{-\frac{1}{2}} = I$}$$

$$= I \tag{By definition of the QR decomposition}$$

We can find the corresponding $B$ by solving $A B = V$, which gives:

$$B = \tilde{Q}^T W^{\frac{1}{2}} V$$

Therefore, this is also a valid solution. The solution above (based on the eigendecomposition; call it $A_{eig}$) and the OP's solution (based on the QR decomposition; call it $A_{QR}$) are related as:

$$A_{eig} = A_{QR} M \quad B_{eig} = M^T B_{QR} \quad \text{where} \quad M = \tilde{Q}^T W^{\frac{1}{2}} V U L^{-\frac{1}{2}}$$

One can show that $M^T M = I$, so this is a unitary transformation (see the note above about infinitely many solutions).

$\endgroup$
  • $\begingroup$ You can do this with modified gram schmidt $\endgroup$ – dave fournier Jan 24 at 5:00
  • $\begingroup$ @davefournier Yes, it seems one could use a version of modified Gram-Schmidt that works in general inner product spaces, rather than the typical Euclidean version. But, the OP asked for a non-recursive solution, which Gram-Schmidt wouldn't satisfy. $\endgroup$ – user20160 Jan 24 at 6:59
  • $\begingroup$ user20160 Not sure why modified gram schmidt is "recursive", but given that why is it more recursive than QR? $\endgroup$ – dave fournier Jan 24 at 16:39
  • $\begingroup$ @user20160 Thank you for your detailed answer. I checked that it indeed provides us orthogonal (as defined in post) matrix. As you mentioned, the form of matrix specified up to rotations. Since I want in the case when all weights are equal to 1 the same result as standard poly() function with QR -- what can I do for it? Could you give a hint how to obtain corresponding rotation? $\endgroup$ – xxxxx Jan 24 at 23:02
  • 1
    $\begingroup$ @xxxxx I took a look again at the solution you proposed based on the QR decomposition. See my edited answer for a proof that it works, and that it's equivalent to the eigendecomposition method (up to rotation & sign flips). So, you should be fine using QR if you want to match the output of poly(). Of course, there's no 'privileged basis' so any unitary transformation is as valid as any other. $\endgroup$ – user20160 Jan 25 at 16:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.