2

I need to solve in Matlab a Least squares problem,

\begin{equation} \begin{bmatrix} x^2 \\ ux \\ u^2 \end{bmatrix}^T\begin{bmatrix} H_{xx} \\ 2H_{xu} \\ H_{uu} \end{bmatrix} = y \end{equation}

Set $z= \begin{bmatrix} x^2 \\ ux \\ u^2 \end{bmatrix}$, $H = \begin{bmatrix} H_{xx} \\ 2H_{xu} \\ H_{uu} \end{bmatrix}$, where $H$ is the vector of unknowns. As I have q=3 unknowns, I need to collect at least $N\geq q$ data samples: \begin{equation} \Phi = [z^1 \ z^2 \ z^3] \end{equation} where the data are collected such that the rank condition is satisfied (rank($\Phi$)=q) and \begin{equation} Y = \begin{bmatrix}y^1 \\ y^2 \\ y^3\end{bmatrix} \end{equation} The Least squares solution is then \begin{equation} H = (\Phi^T\Phi)^{-1}\Phi^T Y \end{equation}

However, the computed H is too big and cannot converge to the true value. This is probably because the matrix $\Phi$ is ill-conditioned (indeed cond(Phi) returns a number of order e+03). What am I missing?

This is the reference I am using: see Algorithm 2

Solved: there was a mistake in the collection of datasample $z^{i}$, so even though the rank condition was satisfied, its inverse was almost singular (indeed the condition number was very high).

Betelgeuse
  • 482
  • So it should be that $z^\top H = y$, so $\Phi^\top H = Y$. So I believe the least squares solution should be $H = (\Phi,\Phi^\top)^{-1} \Phi,Y$ instead. Also what are the singular values of $\Phi$? – Kwin van der Veen May 29 '19 at 12:45
  • for example cond(Phi) =2.8279e+09 and eig(Phi)

    ans =

    -0.5609 + 0.0000i, 0.1212 + 0.2634i, 0.1212 - 0.2634i,

     – Betelgeuse May 29 '19 at 13:18
  • I did the experiments simulating a discrete-time system knowing a vector x of length L and $u = -K x + randn(L,1)$ (to satisfy the rank condition) – Betelgeuse May 29 '19 at 13:23
  • I believe you misinterpreted how to define $z$ (in the paper it is called $\bar{z}_k$). Namely when $x\in\mathbb{R}^2$ and $u\in\mathbb{R}^1$ then according to the paper $$\bar{z}_k=\begin{bmatrix}x_1^2 \ x_1,x_2 \ x_1,u \ x_2^2 \ x_2,u \ u^2\end{bmatrix}.$$ Which would also mean that $H$ has six unknowns. – Kwin van der Veen Jun 01 '19 at 12:08
  • @Kwin van der Veen yes I tried to simplify in the question, assuming just one state. For my problem, I am considering $\bar{z}k$ and H with 6 unknown. Then I am collecting L = 6 data samples of $\bar{z}_k$, $\bar{z}{k+1}$, $r_k$, where the input component $u_k$ has an exploration signal – Betelgeuse Jun 01 '19 at 14:08

1 Answers1

2

Each sample is generated according to $z_i^T H = y_i$. If you have $N$ data points (each one consisting of a three-dimensional vector $z_i$ and an observation $y_i$), you collect them in an $N\times 3$ matrix $$\Phi =\begin{bmatrix} z_1^T\\ \vdots\\z_N^T\end{bmatrix},$$ and an $N\times 1$ vector $$y=\begin{bmatrix} y_1^T\\ \vdots\\y_N^T\end{bmatrix};$$ then, you find the least squares solution and $$\hat H = (\Phi^T \Phi)^{-1}\Phi y.$$

In matlab, you should use the division operation H_hat = Phi\y, which will choose the most appropriate method depending on the matrices (usually it relies on qr factorization).

Riccardo Sven Risuleo
  • 3,056
  • Maybe i collected data in a wrong way. How should I compute the N samples according to the relation $z_i^T H=y_i$ if I don't know the value of $H$? Should I use the true value of H for the simulations? – Betelgeuse May 29 '19 at 15:19
  • I suppose you have excitation variables $z$ and responses $y$ from some data, and you want to find the $H$ that generated them. If you're just testing stuff out, use a true $H$ to generate the data in response to some randomly chosen inputs, then add some random noise to the observations $y$, compute the $\hat H$. As you increase $N$, $\hat H$ will converge to the true $H$. – Riccardo Sven Risuleo May 29 '19 at 15:25
  • Suppose that I have $z_i$ that satisfies the rank condition, I use the true $H$ to generate $y_i$ for N times. If I do then $\hat{H}= pinv(\Phi)Y$, $\hat{H} \neq H$. Same if I generate $z_i$ and $y_i$ ($y$ also depends on $x,u$) based on values $(x_i,u_i)$ (so I don't use $H$ to generate $\Phi,Y$ ), where $u_i = -K x_i+v$, where $v$ random gaussian noise. – Betelgeuse May 29 '19 at 16:04
  • @Betelgeuse But how close is $\hat{H}$ to $H$? Could it just be due to numerical errors? Are you sure that the is a linear relation between $z$ and $y$? – Kwin van der Veen May 29 '19 at 16:36
  • $\hat{H}$ is too far from the true $H$, it has indeed values of order $1.0e+04$. In my problem, considering $z_k$ with elements of $u_k, x_k$ as I defined in my question, $y := x_k^T Q x_k + u_k^TRu_k + H_0 z_{k+1}$ so I have to solve the problem $H_1^T z_k = x_k^T Q x_k + u_k^TRu_k + H_0 z_{k+1}$ ($H_0$ is known, $H_1$ to be found). – Betelgeuse May 29 '19 at 16:43
  • 1
    It could be a persistence of excitation problem; the input data $z_i$ are not informative enough to discriminate all the parameters. In any case, under regularity conditions, least squares is strongly consistent so it converges to the true value w.p. 1. You need to check all the assumptions in your model; I believe the problem is there (for instance, if $z$ depends on $y$ in some way or if you have noise in $z$). – Riccardo Sven Risuleo May 29 '19 at 18:34
  • @Betelgeuse So I assume that $$z_k = \begin{bmatrix}x_k^\top x_k \ u_k^\top x_k \ u_k^\top u_k\end{bmatrix}$$ and is there any relation between $x_k$ and $x_{k+1}$, because currently the term $H_0,z_{k+1}$ doesn't have any correlation? – Kwin van der Veen May 30 '19 at 01:22
  • $x_{k+1}=Ax_k + Bu_{k} + v_k$, where $v_k$ is some random noise. \begin{equation} z_{k+1} = \begin{bmatrix} x_{k+1}^Tx_{k+1} \ u_{k+1}^Tx_{k+1} \ u_{k+1}^Tu_{k+1} end{bmatrix} \end{equation}. – Betelgeuse May 30 '19 at 08:18
  • @Betelgeuse I have done some numerical testing and it seems that $H_1^T z_k = x_k^T Q x_k + u_k^TRu_k + H_0 z_{k+1}$ only seems to hold when $x,u\in\mathbb{R}^2$. – Kwin van der Veen May 30 '19 at 19:54
  • @Kwin van der Veen I am using $x \in \mathbb{R}^2$ and $u \in \mathbb{R}$, I added the reference paper in the question. – Betelgeuse May 31 '19 at 07:34