Feasibility Problem and Similarity Transformations


(Lucas) #1

Hello Everyone,

I’m trying to solve a feasibility problem where all I want to find is a similarity transformation matrix T whose properties should be:

T = T’ (symmetric)

AT == TA’
B == T*C’

where {A,B,C} are known matrices associated with a state-space realization.
Apparently two errors pop-up depending on the parameter values: trivial solution in which case CVX returns NaN for the entries of T or Overdetermined problem and CVX returning the same value. The point is that T is sure to exist in theory, but looks like a numerical solution seems to be elusive so far.
The problem is convex I’m wondering as to what possible trick or equivalence might I be overlooking in especifying this problem.
I’d appreciate some advice from the community.


(Mark L. Stone) #2

You haven’t told us anything about A, B, and C. What is your theoretical basis for being sure a solution exists? Can you show us specific A. B, and C, along with the accompanying theory which shows a solution exists, together with your CVX code and output which fails to find a solution?


(Lucas) #3

Hi Mark, thanks for your reply. Let me elaborate on previous ideas:

Matrix A has entries sitting along the diagonal and some of the superdiagonal and subdiagonal. All entries off this band are zero.
Should look like this:
A = [-a11 0 0 0 0 0
0 -a22 0 0 0 0
0 0 -a33 b 0 0
0 0 -b -a33 0 0
0 0 0 0 -a44 c
0 0 0 0 -c -a44 ]
B = [1
1
2
0
2
0]

C = [c11 c12 c13 c14 c15 c16].

The real A for my problem is sparse and consists of a chain of blocks whose structure equals what I’ve illustrated above whereas the real B has more columns following the pattern given but in some echelon manner (sparse) and C has more rows but is full (not sparse).
The theoretical basis is that every minimal state-space realization (which is the case given above) admits a reciprocal realization thru some similarity transformation. This is a theoretical result from network theory or linear systems theory.

Here goes my code

cvx_begin

variable T(size(A)) symmetric

minimize 0
subject to

T*A ==A’T;
T
B == C’;

cvx_end


(Mark L. Stone) #4

Please show us the output. Can you give a full, reproducible example to work with in which it fails, complete with all input data provided?


(Lucas) #5

Sure! I can provide you with all the data just let me know whether there is a way I can upload the data as a file for the size of these matrices are huge - typing them is not a good idea.


(Mark L. Stone) #6

Can you make a minimal size example which still illustrates your concern?

You ought to be able to easily show the full CVX/solver output.


(Lucas) #7

I’ve tried to to scale it down B

keeping it still meaningful for the applications I have at hand, here they go:


(Lucas) #8

A= [ -1,72E+02 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 ;
0,00E+00 -1,93E+07 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 ;
0,00E+00 0,00E+00 -7,84E+06 1,44E+07 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 ;
0,00E+00 0,00E+00 -1,44E+07 -7,84E+06 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 ;
0,00E+00 0,00E+00 0,00E+00 0,00E+00 -4,57E+06 4,57E+07 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 ;
0,00E+00 0,00E+00 0,00E+00 0,00E+00 -4,57E+07 -4,57E+06 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 ;
0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 -1,72E+02 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 ;
0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 -1,93E+07 0,00E+00 0,00E+00 0,00E+00 0,00E+00 ;
0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 -7,84E+06 1,44E+07 0,00E+00 0,00E+00 ;
0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 -1,44E+07 -7,84E+06 0,00E+00 0,00E+00 ;
0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 -4,57E+06 4,57E+07 ;
0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 0,00E+00 -4,57E+07 -4,57E+06 ]

B = [1 0;
1 0;
2 0;
0 0;
2 0;
0 0;
0 1;
0 1;
0 2;
0 0;
0 2;
0 0];

C= [ 6,94E+00 -2,97E+04 6,88E+04 1,11E+05 1,05E+05 5,28E+03 1,40E+00 -7,40E+01 3,54E+01 1,62E+02 1,61E+02 2,82E+02 ;
1,41E+00 -1,19E+03 -1,66E+02 -2,06E+02 4,49E+02 -4,67E+02 6,72E+00 -3,28E+04 7,65E+04 1,19E+05 1,09E+05 3,84E+03 ];


(Lucas) #10

Sorry about the commas, but I moved from Matlab to Excel and then to this chat.


(Mark L. Stone) #11

When I copy and paste the text, I get

whos
Name Size Bytes Class Attributes

A 12x24 2304 double
B 12x2 192 double
C 2x24 384 double

which is apparently wrong, and dimensionally inconsistent.

And what is the output?


(Lucas) #12

I’ve corrected the commas from excel, now it should be consistent

A= [ -1.72E+02 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 ;
0.00E+00 -1.93E+07 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 ;
0.00E+00 0.00E+00 -7.84E+06 1.44E+07 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 ;
0.00E+00 0.00E+00 -1.44E+07 -7.84E+06 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 ;
0.00E+00 0.00E+00 0.00E+00 0.00E+00 -4.57E+06 4.57E+07 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 ;
0.00E+00 0.00E+00 0.00E+00 0.00E+00 -4.57E+07 -4.57E+06 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 ;
0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 -1.72E+02 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 ;
0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 -1.93E+07 0.00E+00 0.00E+00 0.00E+00 0.00E+00 ;
0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 -7.84E+06 1.44E+07 0.00E+00 0.00E+00 ;
0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 -1.44E+07 -7.84E+06 0.00E+00 0.00E+00 ;
0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 -4.57E+06 4.57E+07 ;
0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 -4.57E+07 -4.57E+06 ];

C= [ 6.94E+00 -2.97E+04 6.88E+04 1.11E+05 1.05E+05 5.28E+03 1.40E+00 -7.40E+01 3.54E+01 1.62E+02 1.61E+02 2.82E+02 ;
1.41E+00 -1.19E+03 -1.66E+02 -2.06E+02 4.49E+02 -4.67E+02 6.72E+00 -3.28E+04 7.65E+04 1.19E+05 1.09E+05 3.84E+03 ];


(Lucas) #13

The dimensionality for B was alright.

For this particular example which is part of a larger system the output is:
Trivial infeasibilities detected; solution determined analytically.
Status: Infeasible
Optimal value (cvx_optval): +Inf


(Mark L. Stone) #14

Yes, it’'s infeasible. It’s not close to feasible.

cvx_begin
variable T(size(A)) symmetric
minimize(norm([T*A-A'*T,T*B-C'],'fro'))
cvx_end

Calling SDPT3 4.0: 139 variables, 60 equality constraints
------------------------------------------------------------

 num. of constraints = 60
 dim. of socp   var  = 117,   num. of socp blk  =  1
 dim. of free   var  = 22
 *** convert ublk to linear blk
********************************************************************************************
   SDPT3: homogeneous self-dual path-following algorithms
********************************************************************************************
 version  predcorr  gam  expon
    NT      1      0.000   1
it pstep dstep pinfeas dinfeas  gap     mean(obj)    cputime    kap   tau    theta
--------------------------------------------------------------------------------------------
 0|0.000|0.000|1.2e+05|7.2e+00|5.4e+04| 2.502306e+03| 0:0:00|5.4e+05|1.0e+00|1.0e+00| chol 1   1
 1|0.011|0.011|1.2e+05|7.2e+00|5.4e+04| 2.499670e+03| 0:0:00|5.4e+05|1.0e+00|9.9e-01| chol 1   1
 2|0.151|0.151|1.0e+05|6.1e+00|4.6e+04| 2.471282e+03| 0:0:00|4.6e+05|1.0e+00|8.4e-01| chol 1   1
 3|0.509|0.509|4.9e+04|3.0e+00|2.3e+04| 2.325334e+03| 0:0:00|2.2e+05|1.0e+00|4.2e-01| chol 1   1
 4|0.967|0.967|1.9e+03|1.2e-01|3.1e+03| 1.675072e+03| 0:0:00|3.4e+03|1.0e+00|1.7e-02| chol 1   1
 5|0.887|0.887|2.2e+02|1.7e-02|3.0e+02| 9.106560e+02| 0:0:00|3.1e+02|1.1e+00|1.9e-03| chol 1   1
 6|0.977|0.977|6.1e+00|4.0e-03|8.5e+00| 1.014258e+03| 0:0:00|1.1e+01|1.1e+00|5.5e-05| chol 1   1
 7|0.988|0.988|1.0e-01|3.4e-03|1.3e-01| 1.016194e+03| 0:0:00|2.1e-01|1.1e+00|9.2e-07| chol 1   1
 8|0.989|0.989|1.7e-03|3.0e-03|1.8e-03| 1.016215e+03| 0:0:00|4.4e-03|1.1e+00|1.5e-08| chol 1   1
 9|0.989|0.989|5.4e-05|1.2e-03|1.6e-04| 1.016215e+03| 0:0:00|8.8e-05|1.1e+00|4.9e-10| chol 1   1
10|0.989|0.989|5.2e-06|5.0e-04|2.3e-05| 1.016215e+03| 0:0:00|2.6e-06|1.1e+00|4.7e-11| chol 1   1
11|0.991|0.991|2.5e-07|2.4e-05|1.1e-06| 1.016215e+03| 0:0:01|1.9e-07|1.1e+00|2.2e-12| chol 1   1
12|0.990|0.990|7.2e-09|6.3e-07|2.9e-08| 1.016215e+03| 0:0:01|1.0e-08|1.1e+00|6.5e-14| chol 1   1
13|0.990|0.990|9.0e-10|1.0e-08|4.9e-10| 1.016215e+03| 0:0:01|3.5e-10|1.1e+00|0.0e+00|
  Stop: max(relative gap,infeasibilities) < 1.49e-08

-------------------------------------------------------------------
 number of iterations   = 13
 primal objective value =  1.01621527e+03
 dual   objective value =  1.01621527e+03
 gap := trace(XZ)       = 4.85e-10
 relative gap           = 4.77e-13
 actual relative gap    = -6.92e-13
 rel. primal infeas     = 8.96e-10
 rel. dual   infeas     = 1.01e-08
 norm(X), norm(y), norm(Z) = 1.4e+03, 7.1e-01, 1.4e+00
 norm(A), norm(b), norm(C) = 4.4e+08, 1.4e+03, 1.0e+00
 Total CPU time (secs)  = 0.53  
 CPU time per iteration = 0.04  
 termination code       =  0
 DIMACS: 9.0e-10  0.0e+00  1.0e-08  0.0e+00  -6.9e-13  2.4e-13
-------------------------------------------------------------------
 
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +1016.22

I think numerically your problem is a little shaky (cond(A) = 2.67 e5), but suspect that is not accounting for the infeasibility.

Anyhow, all we have now is your word that this is supposed to be feasible. This is not really the right forum to adjudicate the correctness of your assertion on theoretical feasibility.


(Mark L. Stone) #15

However, your problem is feasible if you remove the requirement that T be symmetric.

cvx_begin
variable T(size(A))
minimize(norm([T*A-A'*T,T*B-C'],'fro'))
cvx_end


Calling SDPT3 4.0: 145 variables, 0 equality constraints
------------------------------------------------------------

 num. of constraints =  1
 dim. of socp   var  = 145,   num. of socp blk  =  1
 dim. of linear var  =  1
*******************************************************************
   SDPT3: Infeasible path-following algorithms
*******************************************************************
 version  predcorr  gam  expon  scale_data
    NT      1      0.000   1        0    
it pstep dstep pinfeas dinfeas  gap      prim-obj      dual-obj    cputime
-------------------------------------------------------------------
 0|0.000|0.000|4.5e+00|7.4e+00|3.9e+02| 2.408319e+01  0.000000e+00| 0:0:00| chol  1  1 
 1|1.000|1.000|1.5e-06|7.1e-02|3.0e+01| 1.988429e+01 -8.000157e+00| 0:0:00| chol  1  1 
 2|0.996|0.989|5.4e-08|7.8e-03|3.1e-01| 2.167334e-01 -7.730905e-02| 0:0:00| chol  1  1 
 3|0.989|0.989|1.8e-08|7.8e-04|3.4e-03| 2.381520e-03  1.395204e-04| 0:0:00| chol  1  1 
 4|0.989|0.989|7.1e-11|7.9e-05|3.7e-05| 2.617055e-05  1.004342e-04| 0:0:00| chol  1  1 
 5|0.989|0.989|1.6e-10|8.4e-07|4.0e-07| 2.875884e-07  1.075104e-06| 0:0:00| chol  1  1 
 6|0.989|0.989|2.5e-12|9.1e-09|4.5e-09| 3.160068e-09  1.150710e-08| 0:0:00|
  stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
 number of iterations   =  6
 primal objective value =  3.16006808e-09
 dual   objective value =  1.15071028e-08
 gap := trace(XZ)       = 4.48e-09
 relative gap           = 4.48e-09
 actual relative gap    = -8.35e-09
 rel. primal infeas (scaled problem)   = 2.53e-12
 rel. dual     "        "       "      = 9.07e-09
 rel. primal infeas (unscaled problem) = 0.00e+00
 rel. dual     "        "       "      = 0.00e+00
 norm(X), norm(y), norm(Z) = 1.0e+00, 1.2e-08, 1.0e+00
 norm(A), norm(b), norm(C) = 2.0e+00, 2.0e+00, 2.0e+00
 Total CPU time (secs)  = 0.17  
 CPU time per iteration = 0.03  
 termination code       =  0
 DIMACS: 2.5e-12  0.0e+00  9.1e-09  0.0e+00  -8.3e-09  4.5e-09
-------------------------------------------------------------------
 
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +3.16007e-09

Or

cvx_begin
variable T(size(A))
T*A == A'*T;
T*B == C';
cvx_end
Homogeneous problem detected; solution determined analytically.
Status: Solved
Optimal value (cvx_optval): +0

disp(T)

   1.0e+04 *
    0.0007         0         0         0         0         0    0.0001         0         0         0         0         0
         0   -2.9700         0         0         0         0         0   -0.1190         0         0         0         0
         0         0    3.4400    5.5500         0         0         0         0   -0.0083   -0.0103         0         0
         0         0    5.5500   -3.4400         0         0         0         0   -0.0103    0.0083         0         0
         0         0         0         0    5.2500    0.2640         0         0         0         0    0.0225   -0.0233
         0         0         0         0    0.2640   -5.2500         0         0         0         0   -0.0233   -0.0225
    0.0001         0         0         0         0         0    0.0007         0         0         0         0         0
         0   -0.0074         0         0         0         0         0   -3.2800         0         0         0         0
         0         0    0.0018    0.0081         0         0         0         0    3.8250    5.9500         0         0
         0         0    0.0081   -0.0018         0         0         0         0    5.9500   -3.8250         0         0
         0         0         0         0    0.0080    0.0141         0         0         0         0    5.4500    0.1920
         0         0         0         0    0.0141   -0.0080         0         0         0         0    0.1920   -5.4500

(Lucas) #16

Many thanks, Mark. I’m wondering as to why you rewrote the feasibility problem as norm minimization. Whats the rationale?


(Mark L. Stone) #17

I was trying to see how close to feasibility it is.


(Lucas) #18

I supposed you inferred the infeasibility from the narrowing duality gap as well as primal objective value = 1.01621527e+03 dual objective value = 1.01621527e+03 nearly coincident?


(Mark L. Stone) #19

The primal and dual objective values being nearly coincident at 1.01621527e+03 shows that the norm minimization problem was solved, with that as the optimal norm ‘fro’. The fact that it.is not zero shows the infeasibility of the original problem, or more to the point, the “distance” from feasibility.

Running the problem as you originally provided (minimnize 0 statementcan be included or not, same problem) resulted in the CVX/solver claim of infeasibility.