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

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$.

• Why not use vectorization to rewrite as a conventional least-squares problem? – Rodrigo de Azevedo Nov 24 '17 at 8:22
• hi rodrigo, thanks for the suggestion. Can you be more specific by what you mean when you say "vectorization" and "conventional least squares"? – user1237300 Nov 24 '17 at 15:01
• Vectorization. Conventional least squares is to minimize $\| A x - b\|_2^2$, which produces the normal equations. – Rodrigo de Azevedo Nov 24 '17 at 15:09

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)$.

• Thank you so much for the detailed response, Brian! It did not occur to me that you can rewrite the unknown C into a vector. This is why I didn't know how to proceed with regularization of a matrix. As a side question, I read about Matlab's pinv which is supposed to give the minimum 2 norm solution. Hence, in $C=Y∖(Q/Z)$ above, I replaced all slashes with pinv and putting the matrix transposes as necessary. In other words, I computed the minimum norm solution of solving the least squares problem $Q/Z$ then using this minimum norm solution, I obtained the minimum norm solution of $Y\(Q/Z)$. – user1237300 Nov 26 '17 at 20:23
• I know that the above is just a heuristic, but does this 2 step process of obtaining least norm solutions make sense? – user1237300 Nov 26 '17 at 20:24
• It is not clear to me that the two-step process gives a least-squares solution or that if it does give a least-squares solution that it also gives a minimum norm least squares solution. – Brian Borchers Nov 26 '17 at 20:28