Solving for $C$ in $Q = YCZ$ using least squares in Matlab

Question

I am trying to solve for the matrix $C$ in $Q = YCZ$ in matlab. I have preliminary results but they don't seem realistic. Here, $Q$ is $n \times m-1$, $Y$ is $n \times p$, $C$ is $p \times m$ and $Z$ is $m \times m-1$.

The value of $n$ is very large, say $n = 10000$. On the other hand, $m$ and $p$ are of the same order of magnitude. In the case I considered, $m=201, p =175$.

Among all matrices above, $Z$ is the only one with a seemingly unique structure. By eliminating the first row of $Z$, it becomes upper triangular and invertible.

In matlab, I tried to solve for $C$ using: $C = Y\backslash (Q/Z)$. However, I ran into the following issues:

1) The last column of $C$ is all zero which is unrealistic in my problem. It turns out that the last column of $Q/Z$ is also a zero vector which is probably the cause of this. Are there suggestions how to compute $Q/Z$ without having the last column being 0?

2) If I want to find $C$ such that every row of $C$ has minimum norm, how can I modify my problem in matlab? Based on the motivation below, every row of $C$ corresponds to values of a function on the unit interval. It is highly preferable for these functions to not be highly oscillatory in contrast to what I'm seeing now.

Thanks for any suggestions on how to proceed.

The motivation of the problem is as follows: let $Q(x)$ be a stochastic process defined as $Q(x) = \int_0^x P(s)\,ds$ where $P(s)$ is another stochastic process. Suppose that $P(s) = q_1(s)Y_1 + \dots + q_n(s) Y_n$ where $Y_1,\dots,Y_n$ are random variables and $q_i(s)$ are deterministic functions of $s$. It will be assumed that the domain of $P(s)$ and $Q(x)$ is the interval $[0,1]$. Suppose that I have observations of $Q(x)$ at discrete time points and that this information is stored in $Q$. Each row of this matrix is an observation while every column corresponds to a point in the spatial discretization. In addition, suppose that I have observations of $Y_1,\dots,Y_n$ stored in a matrix $Y$ where every row is an observation and every column corresponds to $Y_i$. My interest is to estimate the deterministic functions $q_1(s),\dots,q_n(s)$. To do this, I use a Finite element basis, i.e. example, such that $q_i(s) = \sum_{k=1}^m c_i^k \phi_k (s)$ where $\phi_k$ are the basis functions and $c_i^k$ are the nodal values. The matrix $C$ is then obtained such that $C_{ik} = c_i^k$ using notation above and is the quantity of interest. The definition of $Z$ follows from the relationship $Q(x) = \int_0^x P(s) \,ds$.

Answer 1

You can reformulate this problem as a conventional linear least squares problem as follows:

First, write $YCZ$ as

$YCZ=\sum_{i=1}^{p} \sum_{j=1}^{m} C_{i,j} (Y_{:,i}Z_{j,:})$

Next, define the $\mbox{vec}$ operator, which stacks column $1$ of a matrix over columns $2$, $3$, $\ldots$ to form a vector. We'll adopt the convention that lowercase letters are used to denote the vector version of the corresponding matrix which is denoted by an uppercase letter.

Now, we can write

$\mbox{vec}(YCZ)=Ac$

where $c=\mbox{vec}(C)$ and $A$ is

$A=\left[ \mbox{vec}(Y_{:,1}Z_{1,:}) \; \mbox{vec}(Y_{:,1}Z_{2,:}) \; \ldots \; \mbox{vec}(Y_{:,p} Z_{m,:}) \right]$

Your least squares problem is now

$\min \| Ac-q \|_{2}^{2}$.

Although this problem appears to be overdetermined because $A$ has $n(m-1)$ rows and only $pm$ columns, it's actually column rank deficient because of the special structure of $A$. This means that there will be infinitely many least squares solutions. MATLAB's \ operator gives a least squares solution with c=A\q, but provides no control over which of the infinitely many least squares solutions you'll get.

Rather, you'll want to include some kind of regularization term. For example, if you want to minimize the Frobenius norm of $C$, you could minimize

$\min \| Ac-q \|_{2}^{2} +\lambda^{2} \| c \|_{2}^{2}$

where $\lambda$ is a small regularization parameter.

This can be written as

$\min \| \bar{A}c - \bar{q} \|_{2}^{2}$

where

$\bar{A}=\left[ \begin{array}{c} A \\ \lambda I \end{array} \right]$

and

$\bar{q}=\left[ \begin{array}{c} q \\ 0 \end{array} \right]$.

Other kinds of quadratic regularization are relatively easy to obtain by replacing $\lambda I$ with an appropriate matrix.

In actually solving this linear least squares problem, you'll want to avoid forming the very large (and presumably dense) matrix $A$. You can use an iterative method such as LSQR (included in MATLAB) for the least squares problem and take advantage of the special structure of $A$ to provide your own routines for multiplying $A$ times a vector and $A^{T}$ times a vector. For example, to multiply $Au$, simply take $\mbox{vec}(YUZ)$ where $u=\mbox{vec}(U)$.