Message Boards Message Boards

Solve f(X)=0 where X is a matrix and f(X) is a complex matrix?

GROUPS:

Why is the solution a {2,2,2} dimension instead of {2,2}? (Thank to Bill Simpson)

Clear[a, b, c, d];
X = {{a, b}, {c, d}};
h1 = RandomComplex[{2 + I, 10 + 20 I}, {2, 2}];(*MIMO channel*)
h2 = RandomComplex[{2 + I, 10 + 20 I}, {2, 2}];
myu = RandomReal[1, 2]; (*probabilities of the realization of the MIMO channels*)
Lambda = 0.3;
FX1 = myu[[1]]*ConjugateTranspose[h1].Inverse[
     ConjugateTranspose[h1.X.h1]/0.2 + IdentityMatrix[2]].h1;
FX2 = myu[[2]]*ConjugateTranspose[h2].Inverse[
     ConjugateTranspose[h2.X.h2]/0.2 + IdentityMatrix[2]].h2;
FXlambda = FX1 + FX2;
Map[X /. # &, 
 NSolve[myu[[1]]*
     ConjugateTranspose[h1].Inverse[
       ConjugateTranspose[h1.X.h1]/0.2 + IdentityMatrix[2]].h1 - 
    Lambda*IdentityMatrix[2] == RandomComplex[{0, 0}, {2, 2}], {a, b, c, d}]]
POSTED BY: Massa Ndong
Answer
27 days ago

(1) Please use the formatting tools above the pane to put code into a proper code block. The first icon, with the <>, is used to make a code block.

(2) Use correct Wolfram Language comments so that they need not be removed in a cut-and-paste in order to actually run the code.

(3) There are two solutions and each is a 2x2 matrix, all in accordance with the dimensionality and degree of the system.

POSTED BY: Daniel Lichtblau
Answer
27 days ago

Hi Daniel,

Thanks for the comments. How to better solve my f(X)=O with an initial guess? I'm looking for real solutions on the diagonal of X.

Clear[a, b, c, d];
X = {{a, b}, {c, d}};
h1 = RandomComplex[{2 + I, 10 + 20 I}, {2, 2}];(*MIMO channel*)
h2 = RandomComplex[{2 + I, 10 + 20 I}, {2, 2}];
myu = RandomReal[1, 
   2];(*probabilities of the realization of the MIMO channels*)
Lambda = 0.3;
FX1 = myu[[1]]*
   ConjugateTranspose[h1].Inverse[
     X.h1.ConjugateTranspose[h1]/0.2 + IdentityMatrix[2]].h1;
FX2 = myu[[2]]*
   ConjugateTranspose[h2].Inverse[
     X.h2.ConjugateTranspose[h2]/0.2 + IdentityMatrix[2]].h2;
FXlambda = FX1 + FX2;
(*Map[X/.#&,NSolve[FX1+FX2\[Equal]RandomComplex[{0,0},{2,2}],{a,b,c,d}\
]]*)
FindRoot[FXlambda == 
  RandomComplex[{0, 0}, {2, 2}], {{a, 0.01}, {b, 0.12}, {c, 0.04}, {d,
    0.14}}]
POSTED BY: Massa Ndong
Answer
20 days ago

This

Clear[a, b, c, d];
X = {{a, b}, {c, d}};
h1 = RandomComplex[{2 + I, 10 + 20 I}, {2, 2}];(*MIMO channel*)
h2 = RandomComplex[{2 + I, 10 + 20 I}, {2, 2}];
myu = RandomReal[1, 2];(*probabilities of the realization of the MIMO channels*)
Lambda = 0.3;
FX1 = myu[[1]]*ConjugateTranspose[h1].Inverse[X.h1.ConjugateTranspose[h1]/0.2 + IdentityMatrix[2]].h1;
FX2 = myu[[2]]*ConjugateTranspose[h2].Inverse[X.h2.ConjugateTranspose[h2]/0.2 + IdentityMatrix[2]].h2;
FXlambda = FX1 + FX2;
min = FindMinimum[Simplify[Norm[Diagonal[FXlambda - RandomComplex[{0, 0}, {2, 2}]]]],
   {{a, 0.01}, {b, 0.12}, {c, 0.04}, {d, 0.14}}]

rapidly looks for approximate real solutions on the diagonal given your initial guess and returned

{5.27366*10^-16, {a -> 9.82904*10^14, b -> 5.07643*10^15, c -> 9.80715*10^14, d -> 3.55582*10^15}}

FXlambda /. min[[2]] is then

{{-8.892*10^-17 - 1.30743*10^-16 I, 1.36888*10^-16 - 7.85417*10^-17 I},
 {3.14381*10^-16 + 5.59592*10^-16 I, -4.8582*10^-16 + 1.30743*10^-16 I}}

If you need more precision and you can justify this then you may want to make the following changes

Clear[a, b, c, d];
X = {{a, b}, {c, d}};
h1 = RandomComplex[{2 + I, 10 + 20 I}, {2, 2}, WorkingPrecision -> 32];(*MIMO channel*)
h2 = RandomComplex[{2 + I, 10 + 20 I}, {2, 2}, WorkingPrecision -> 32];
myu = RandomReal[1, 2, WorkingPrecision -> 32];(*probabilities of the realization of the MIMO channels*)
Lambda = 3/10; (* 0.3 has only 16 digits of precision, use rationals or more precise decimal constants *)
FX1 = myu[[1]]*ConjugateTranspose[h1].Inverse[X.h1.ConjugateTranspose[h1]/(2/10) + IdentityMatrix[2]].h1;
FX2 = myu[[2]]*ConjugateTranspose[h2].Inverse[X.h2.ConjugateTranspose[h2]/(2/10) + IdentityMatrix[2]].h2;
FXlambda = FX1 + FX2;
min = FindMinimum[Simplify[Norm[Diagonal[FXlambda - RandomComplex[{0, 0}, {2, 2}]]]],
   {{a, 1/100}, {b, 12/100}, {c, 4/100}, {d, 14/100}}, WorkingPrecision -> 32, AccuracyGoal -> 32]

which then returned

{9.6593382693239849291748344551116*10^-59,
 {a -> 2.7682397312684356156723753595322*10^57, 
  b -> 7.2319582199103338316464728892929*10^57, 
  c -> 2.0546745021764832416496939588821*10^57, 
  d -> 2.9972363123811069167375347579589*10^58}}

I think these results give hints about the behavior of your system.

As always, you may need to choose your initial guesses near the solutions you are looking for if you are hoping to find particular solutions. You can also try many different initial guesses to see if this can locate more than one solution.

POSTED BY: Bill Simpson
Answer
19 days ago

Hi Bill, I do appreciate your comments. I'm looking for something where c=b. I'm working on a loop structure to build the vector FX containing the function FXi, i_1,2...N and adjust Lamda to keep up with the constraint a+d=Pmax. I did the following but I'm getting a kind of comment.

Clear[a, b, c, d];
X = {{a, b}, {c, d}};
Nt = 2; (*Number of transmit antennas*)
Nr = 2; (*Number of receive antennas*)
N0 = 2/10; (*Noise energy*)
         (*MIMO channel*)
h1 = RandomComplex[{45 + I, 10 + 20 I}, {Nr, Nt}, 
  WorkingPrecision -> 42];
h2 = RandomComplex[{105 + I, 100 + 20 I}, {Nr, Nt}, 
   WorkingPrecision -> 42];
h3 = RandomComplex[{20/10 + 1/100 I, 1 + 2 I}, {Nr, Nt}, 
   WorkingPrecision -> 42];
         (*Singular values*)
SV1 = SingularValueList[h1, Min[Nr, Nt]];
SV2 = SingularValueList[h2, Min[Nr, Nt]];
SV3 = SingularValueList[h3, Min[Nr, Nt]];
                       (*Condition numbers*)
ConNum1 = 10*Log[10, Max[SV1]/Min[SV1]];
ConNum2 = 10*Log[10, Max[SV2]/Min[SV2]];
ConNum3 = 10*Log[10, Max[SV3]/Min[SV3]];

(*probabilities of the realization of the MIMO channels*)
myu = RandomReal[1, 3, WorkingPrecision -> 42];
(*Lagrange multiplier*)
Lambda = 3/10;
(* (FXlambda-lambda*I) is the Lagrangian, where I is the identity \
matrix of dimension Nt *)
FX1 = myu[[1]]*
   ConjugateTranspose[h1].Inverse[
     X.h1.ConjugateTranspose[h1]/N0 + IdentityMatrix[Nr]].h1;
FX2 = myu[[2]]*
   ConjugateTranspose[h2].Inverse[
     X.h2.ConjugateTranspose[h2]/N0 + IdentityMatrix[Nr]].h2;
FX3 = myu[[3]]*
   ConjugateTranspose[h3].Inverse[
     X.h3.ConjugateTranspose[h3]/N0 + IdentityMatrix[Nr]].h3;
FXlambda = FX1 + FX2 + FX3 - Lambda*IdentityMatrix[Nt];

(*Inititial condition is X/Nt*)

min = FindMinimum[
  Simplify[Norm[
    Diagonal[FXlambda - RandomComplex[{0, 0}, {2, 2}]]]], {{a, 
    1/Nt}, {b, 0}, {c, 0}, {d, 1/Nt}}, WorkingPrecision -> 42, 
  AccuracyGoal -> 42]

Which returns

FindMinimum::precw: The precision of the argument function (\[Sqrt](Abs[((1.45407*10^-13+0. I)+(3.35538*10^-9-3.48721*10^-12 I) c+<<15>>+Power[<<2>>] Plus[<<6>>])/(((1.41108*10^-8+0. I)+Times[<<2>>]+Times[<<2>>]+Times[<<2>>]+Times[<<3>>]+Times[<<2>>]+Times[<<3>>]) (<<1>>) (<<1>>))]^2+Abs[((1.67902*10^-13+0. I)+<<17>>+Power[<<2>>] Plus[<<6>>])/(((1.41108*10^-8+0. I)+<<5>>+Times[<<3>>]) (<<1>>) (<<1>>))]^2)) is less than WorkingPrecision (42.`).

FindMinimum::cvmit: Failed to converge to the requested accuracy or precision within 100 iterations.

{7.73626635300012230044601812134030093708021*10^-21, {a -> 
1.14497507023363022694010787328331353072246, 
b -> -0.0187196526797612369572634135990579479468498, 
c -> -0.0187196526797612175981886823921262272387266, 
d -> 1.15134717083103909872642135639044335920618}}

The solution looks good. If I could reach a higher accuracy while maintaining the b=c, that would be great.

POSTED BY: Massa Ndong
Answer
19 days ago

Hi Bill,

Why do you need the 'Diagonal' in

min = FindMinimum[Simplify[Norm[Diagonal[FXlambda - RandomComplex[{0, 0}, {2, 2}]]]],
   {{a, 1/100}, {b, 12/100}, {c, 4/100}, {d, 14/100}}, WorkingPrecision -> 32, AccuracyGoal -> 32]
POSTED BY: Massa Ndong
Answer
13 days ago

I included 'Diagonal' because I was attempting to interpret your 'I'm looking for real solutions on the diagonal of X.' as your primary goal.

If that was not your goal then perhaps you could try to explain more clearly what you are looking for.

POSTED BY: Bill Simpson
Answer
12 days ago

I'm looking for an X with real elements on its diagonal where X is the solution of the equation:

FXlambda==RandomComplex[{0, 0}, {2, 2}]
POSTED BY: Massa Ndong
Answer
12 days ago

This

RandomComplex[{0, 0}, {2, 2}, WorkingPrecision->64]

is just {{0,0},{0,0}} because you are looking for a matrix with each value between 0 and 0 and thus can be discarded.

This

Clear[a, b, c, d];
X = {{a, b}, {c, d}};
h1 = RandomComplex[{2 + I, 10 + 20 I}, {2, 2}, WorkingPrecision->64];(*MIMO channel*)
h2 = RandomComplex[{2 + I, 10 + 20 I}, {2, 2}, WorkingPrecision->64];
myu = RandomReal[1, 2, WorkingPrecision->64];(*probabilities of the realization of the MIMO channels*)
Lambda = 3/10;(*0.3 has only 16 digits of precision,use rationals or more precise decimal constants*)
FX1 = myu[[1]]*ConjugateTranspose[h1].Inverse[X.h1.ConjugateTranspose[h1]/(2/10)+IdentityMatrix[2]].h1;
FX2 = myu[[2]]*ConjugateTranspose[h2].Inverse[X.h2.ConjugateTranspose[h2]/(2/10)+IdentityMatrix[2]].h2;
FXlambda = FX1 + FX2

gives you a very very large random complex matrix over a,b,c,d. I could show you the particular random matrix I generated from this particular evaluation, which would take many screens here, but you can just as easily generate your own random complex matrix.

This

min = FindMinimum[Simplify[Norm[Diagonal[FXlambda]]], {{a,1/100}, {b,12/100}, {c,4/100}, {d,14/100}},
   WorkingPrecision->64, AccuracyGoal->64, MaxIterations->10^4]

returns at this particular moment with these particular random values

{1.409337889826012221043454568743455506805196535480069698590240944*10^-91,
{a -> -3.616219280166386544963592969780528069042039293961884040397257279*10^90,
b -> -1.047395840122292892172381403380611274109197612363848011420483615*10^91,
c -> -4.265655941117427663894689038211955208567032698493186994406768870*10^90, 
d -> -1.933152221127828747959356370666759451780057753648769872401242173*10^91}}

with occasional complaints about convergence and accuracy which you can work on overcoming.

This

X /. min[[2]]

returns this

{{-3.616219280166386544963592969780528069042039293961884040397257279*10^90,
   -1.047395840122292892172381403380611274109197612363848011420483615*10^91},
 {-4.265655941117427663894689038211955208567032698493186994406768870*10^90,
   -1.933152221127828747959356370666759451780057753648769872401242173*10^91}}

which is your desired matrix X with Real elements on the diagonal.

This

FXlambda /. min[[2]]

returns

{{-1.228521870908407132391231192817176855185507203455432182381434*10^-91 +
 1.97262984247407592288931878416167615378333567798851071480550*10^-92 I, 
 4.44888691826143307082987838124192490988725917094492380233716*10^-92 +
 2.33797978252268798471808096323191424766051989096546366847254*10^-92 I},
 {9.47225280258610185565133151101475581141622957605503796821495*10^-92 -
 3.01547536822613285530732469452419099256699443356422071748916*10^-92 I,
-6.31776804330895050336828430370877474291208836448824238836773*10^-92
- 1.97262984247407592288931878416167615378333567798851071480550*10^-92 I}}

which demonstrates that the solution X gives a result FXlambda=={{0,0},{0,0}}, at least within the tiny numerical error associated with your 64 digit random complex original matricies.

Or if you don't want to see the fine details of the solution then

Chop[FXlambda /. min[[2]]] == RandomComplex[{0, 0}, {2, 2}, WorkingPrecision -> 64]

will return

True

I believe this is one solution to your problem. Changing the initial values for a,b,c,d in the FindMinimum may or may not find different solutions.

If you want c==b then

Clear[a, b, d];
X = {{a, b}, {b, d}};
h1 = RandomComplex[{2 + I, 10 + 20 I}, {2, 2}, WorkingPrecision -> 64];(*MIMO channel*)
h2 = RandomComplex[{2 + I, 10 + 20 I}, {2, 2}, WorkingPrecision -> 64];
myu = RandomReal[1, 2, WorkingPrecision->64];(*probabilities of the realization of the MIMO channels*)
Lambda = 3/10;(*0.3 has only 16 digits of precision,use rationals or more precise decimal constants*)
FX1 = myu[[1]]*ConjugateTranspose[h1].Inverse[X.h1.ConjugateTranspose[h1]/(2/10)+IdentityMatrix[2]].h1;
FX2 = myu[[2]]*ConjugateTranspose[h2].Inverse[X.h2.ConjugateTranspose[h2]/(2/10)+IdentityMatrix[2]].h2;
FXlambda = FX1 + FX2;
min = FindMinimum[Simplify[Norm[Diagonal[FXlambda]]], {{a, 1/100}, {b, 12/100}, {d, 14/100}},
   WorkingPrecision -> 64, AccuracyGoal -> 64, MaxIterations -> 10^4]
X /. min[[2]]

returns

{{1.042857763942349304591968267547967620603583294015470931877677227*10^164,
  5.143214731304323107915050538300431747411562194323159341662060110*10^164},
  {5.143214731304323107915050538300431747411562194323159341662060110*10^164, 
  1.470119158978919648478190964930733450490851469062127614340916143*10^164}}

and

Chop[FXlambda /. min[[2]]] == RandomComplex[{0, 0}, {2, 2}, WorkingPrecision -> 64]

returns

True

If I make the following changes

Clear[a, b, d];
X = {{a, b}, {b, d}};
Nt = 2;(*Number of transmit antennas*)
Nr = 2;(*Number of receive antennas*)
N0 = 2/10;(*Noise energy*)(*MIMO channel*)
h1 = RandomComplex[{45 + I, 10 + 20 I}, {Nr, Nt}, WorkingPrecision -> 128];
h2 = RandomComplex[{105 + I, 100 + 20 I}, {Nr, Nt}, WorkingPrecision -> 128];
h3 = RandomComplex[{20/10 + 1/100 I, 1 + 2 I}, {Nr, Nt}, WorkingPrecision -> 128];
(*Singular values*)
SV1 = SingularValueList[h1, Min[Nr, Nt]];
SV2 = SingularValueList[h2, Min[Nr, Nt]];
SV3 = SingularValueList[h3, Min[Nr, Nt]];
(*Condition numbers*)
ConNum1 = 10*Log[10, Max[SV1]/Min[SV1]];
ConNum2 = 10*Log[10, Max[SV2]/Min[SV2]];
ConNum3 = 10*Log[10, Max[SV3]/Min[SV3]];
(*probabilities of the realization of the MIMO channels*)
myu = RandomReal[1, 3, WorkingPrecision -> 128];
(*Lagrange multiplier*)
Lambda = 3/10;
(*(FXlambda-lambda*I) is the Lagrangian,where I is the identity matrix of dimension Nt*)
FX1 = myu[[1]]*ConjugateTranspose[h1].Inverse[X.h1.ConjugateTranspose[h1]/N0+IdentityMatrix[Nr]].h1;
FX2 = myu[[2]]*ConjugateTranspose[h2].Inverse[X.h2.ConjugateTranspose[h2]/N0+IdentityMatrix[Nr]].h2;
FX3 = myu[[3]]*ConjugateTranspose[h3].Inverse[X.h3.ConjugateTranspose[h3]/N0+IdentityMatrix[Nr]].h3;
FXlambda = FX1 + FX2 + FX3 - Lambda*IdentityMatrix[Nt];
(*Inititial condition is X/Nt*)
min = FindMinimum[Simplify[Norm[Diagonal[FXlambda]]], {{a, 1/Nt}, {b, 0}, {d, 1/Nt}}, 
  WorkingPrecision -> 128, AccuracyGoal -> 64, MaxIterations -> 10^4]

then I get additional accuracy, but still warnings that you should explore. And you should inspect

FXlambda /. min[[2]]

carefully after you have modified your problem in this way because it appears that the result is no longer approximately zero.

You might think carefully, exactly what it is that you really want to find the minimum of in your FXlambda. Is it the Norm of the Diagonal? Is it the Total of the Norm of every element in FXlambda?

POSTED BY: Bill Simpson
Answer
12 days ago

Thank for the details. I'm trying to solve FXlambda==0 where

FXlambda = FX1 + FX2 + FX3 - Lambda*IdentityMatrix[Nt];

My understanding is that is equivalent to solving a system of equations made of the elements of the matrix FXlambda. Thus I think minimizing the total of the Norm of every element in FXlambda would be what I want.

POSTED BY: Massa Ndong
Answer
11 days ago

Then your last example

min = FindMinimum[Simplify[Norm[Flatten[FXlambda]]], {{a, 1/Nt}, {b, 0}, {d, 1/Nt}}, 
  WorkingPrecision -> 128, AccuracyGoal -> 64, MaxIterations -> 10^4];
FXlambda /. min[[2]]

or

min = FindMinimum[Simplify[Norm[FXlambda[[1, 1]]] + Norm[FXlambda[[1, 2]]] + 
 Norm[FXlambda[[2, 1]]] + Norm[FXlambda[[2, 2]]]], {{a, 1/Nt}, {b, 0}, {d, 1/Nt}}, 
  WorkingPrecision -> 128, AccuracyGoal -> 64, MaxIterations -> 10^4];
FXlambda /. min[[2]]

do not seem to find a minimum that is almost exactly zero, but perhaps this is close enough.

Or perhaps that starting point led to a local minimum.

POSTED BY: Bill Simpson
Answer
11 days ago

Group Abstract Group Abstract