Numerical solution to complicated implicit matrix equation

Hi,

I would like to ask if someone can give me some advice on how to solve (numerically) a complicated matrix equation system which appeared when trying to solve a constrained optimization problem.

The equation system is as follows:

Defining: A = H*inv(G'*G + L)*H' and C = inv(A+B^2)*A*inv(A+B^2)

       S = inv(H'*inv(B)^2*H + G'*G)*H'*inv(B)

G and H are two constant complex-valued matrices of size N*r-by-M*t. L and B are diagonal matrices, and their diagonal entries are precisely the independent variables of the problem. Their sizes are M*t-by-M*t and N*r-by-N*r respectively. However, L containts only M independent variables (each variable appears t times in the diagonal). Mathematically, L = kron(diag(lambda),eye(t)) where lambda is a vector of length M with the independent variables of L. The independent variables are all real-valued.

I want to find the values of those independent variables such that:

C(i,i) == alpha(i)/(B(i,i)^4) for all diagonal entries of C and B (index i) kron(eye(M),ones(1,t))*abs(S).^2 == P*ones(M,1)

Both alpha and P are also constants in the problem. This leaves me with a system of M+N*r non-linear equations with M+N*r variables. As an extra piece of information (in case it is helpful), the M independent variables within matrix L are Lagrange multipliers coming from the constrained optimization problem. The solution to that equation exists and it is unique, since it comes from a convex optimization problem. However, the matrix sizes tend to be big (around 300-by-300) and together with the fact that there are many matrix inverses in the formulation, the problem exhibits a pathological behavior.

All the "typical" functions such a fsolve have been unable to converge and I am completely out of ideas for now. I know the problem is pretty complicated but if any of you know of something worth trying, a function or even a numerical algorithm able to solve (or approximate the solution), I would be really grateful.

Answer 1

[Star Strider]

First, I suggest you recode these lines:

         A = H*inv(G'*G + L)*H' 

         C = inv(A+B^2)*A*inv(A+B^2)

         S = inv(H'*inv(B)^2*H + G'*G)*H'*inv(B)

as:

         A = H/(G'*G + L)*H' 

         C = (A+B^2)\A/(A+B^2)

         S = (H'/B^2\H + G'*G)\H'/B

Thhe ‘/’ ‘\’ operators eliminate the necessity of calculating the explicit inverses. (I believe I recoded those correctly. Be sure to check them to be sure.)

[Daloal]

Thanks for your comments.

Actually, I should have clarified in my previous post that the "code" I wrote was just to show the shape of the equations ( I actually wanted to use latex but this is my first post here and I don't know how to insert tex). In the implementation I have tried, the equations are implemented as you suggested, without calling the inv function. When the matrix sizes are small (below 10 approximately), avoiding the inv function leads to proper convergence.

Sadly, when the sizes grow an order of magnitude, as it is required in my application, the method breaksdown even when programmed like this.

[Star Strider]

I'm not quite sure what the structures of your matrices are, but you might consider the sparse matrix functions: http://www.mathworks.com/help/techdoc/ref/sparse.html.

I haven't needed to use them in my applications so I have no experience with them, but they might help you in yours.

[Daloal]

The structure is pretty bad in that sense too. They are not sparse in nature (except in the two diagonal matrices with the independent variables, the others have all its entries non-zero).

Moreover, the entries are all of them complex numbers with a very small modulus (in the order of 1e-6 to 1e-10) since they represent signal propagation in an urban environment (lot of losses). When computing the inverses, this gives a lot of trouble too in terms of convergence of numerical solvers.

[Star Strider]

It might be possible to scale your data and then rescale the results of your calculations later. I have had to do this from time to time (with physiolgical data when I was dealing with orders-of-magnitude differences in different sets of data) but not with complex data and not in your applications.

Other than that, all I can do is give you my sympathies.