First kind of fredholm equation Ax=b than


(shine_wade) #1

First kind of fredholm equation Ax=b then perform FFT on the time domain signal,change the equation to AX=B, A&B is in complex domain ,is it feasible to use cvx, or please can give suggestion.


(Mark L. Stone) #2

You need to state more clearly the optimization or feasibility problem you wish to solve. Is it convex? Please read Why isn't CVX accepting my model? READ THIS FIRST!


(shine_wade) #3

Thank you for your immediate reply,it is ill-posed problem,Regularization(like tikhonov is applicable image ).Can be simplified to A(t,z)*x(z)=b(t),where t is time,z is position ;I use FFT for this equation ,in order to solve x(z), to use some tools package would be easy for me


(shine_wade) #4

after FFT(it’s usefule for my work) i want to solve x(z) in complex domain@Mark_L_Stone,:grinning:


(Mark L. Stone) #5

I still am not clear on what you want to do.

If xis a real or complex vector and is the only decision (CVX) variable; and b, \lambda, and x^* are all numerical input data, then the objective function can be entered into CVX, using square_pos(norm(...)) . If there are only linear equality constraints on x, even if they are complex variable constraints, they can be straightforwardly entered into CVX, and the problem can be solved use any of CVX;s available solvers. Please read the CVX Users’ Guide http://cvxr.com/cvx/doc/ to learn the details.


(shine_wade) #6

thankyou so much ,and i still need your help ,error as follow :
% A is complex matrics (mn); x is a real vector(n1) need to gain;b is a
%complex vector (m*1)

cvx_begin
variables x(n,1) A(m,n) b(m,1) complex
expression Objective
Objective = square_pos(norm((A*x(:,1)-b(:,1)),2) + lambda^2 * norm(x(:,1),2)); 
minimize(sum(Objective))   
cvx_end

error:错误使用 * (line 109)
Disciplined convex programming error:
Invalid quadratic form: not a scalar.

error in : Objective = square_pos(norm((A*x(:,1)-b(:,1)),2) + lambda^2 * norm(x(:,1),2));


(Mark L. Stone) #7

I don;t know what problem you ran. The error message does not match the code you posted . PHI appears in the error message, but not in the code you posted.

If possible, please show a compete, preferably small,reproducible example (complete with all input data,ready for readers to copy and paste your code into their own MATLAB session) which exhibits the error you are receiving.

Also note that you didn’t square (using square_pos) ,norm(x(:,1),2)), but the mathematical display in a previous post shows that term being squared.


(shine_wade) #8
tic;
close all
format short
thmcd=0.22; 	
dens=910;   	
c=1900;     		
A=1.54*10^-4;  	    
yita=0.1;   		
q=0.005;    		
tcth=0.2;   		
D=1.2724*1e-7;
e1=10*1e6;
d=9.8*1e-6;      	
multi=1;
D=D*multi;   				
StartPosition=9.3e-8*exp(0.8e5*d)+1.2e-6;

tc=(d)^2/(pi^2*D);  			
T0=yita*q/(c*dens*A*d);   
az=1.35*10^-4;  		
aep=2.45*az;   		
ep0=8.8542*10^-12;  	
epr=2.2;    			
ts=1e-7;       		
tstart=1e-7;   		
tstop=0.005;   		
t=tstart:ts:tstop;
z=linspace(0,d,100); 

fmin=1/tstop;
fmax=4*10^5;               	
E1=linspace(0,10,25);
E2=linspace(10,10,length(z)-56);
E3=linspace(10,0,31);
E=[E1 E2 E3];
G=(aep-az)*ep0*epr*E;
n=1:100;
T=T0*ones(length(z),1)*exp(-t/tcth)...             	 	
    +2*T0*cos(pi*z'*n/d)*exp(-(n.^2)'*t/tc);
dTdt=-T0/tcth*ones(length(z),1)*exp(-t/tcth)...     		
    +2*T0*(ones(length(z),1)*(-n.^2/tc)).*cos(pi*z'*n/d)*exp(-(n.^2)'*t/tc);
Gtemp=zeros(1,length(z)-1);
Tt1=zeros(length(z)-1,1);
dTdttemp=zeros(length(z)-1,length(t));
for k=1:length(z)-1
    Gtemp(k)=(G(k)+G(k+1))/2;
    Tt1(k)=(T(k,1)+T(k+1,1))/2;
    dTdttemp(k,:)=(dTdt(k,:)+dTdt(k+1,:))/2;
end

If=A/(length(z)-1)*Gtemp*dTdttemp;
Qt1f=A/(length(z)-1)*Gtemp*Tt1;
f=fmin:fmin:fmax;
nf=length(f);
Yf=ts*If*exp(-1i*2*pi*t'*f)+Qt1f; 
Yfr=real(Yf);
Yfi=imag(Yf);
w=2*pi*f;
k=(1+1i)*(sqrt(w./(2*D)))';
ks=k./sinh(k.*d);
Tw0=yita*q/(c*dens*A);  
Tw=Tw0*zeros(size(f,2),size(z,2))+ks.*cosh(k.*(d-z));

%Kx=b
m1=size(K,1);  
n1=size(K,2);  
K=Tw.*(aep-az)*ep0*epr;
b=Yf';
%K;   %m1 x n1 matrix    
lambda = 1.045929285995968e-09;

%Use cvx
tic
cvx_begin
variables x(n1,1) K(m1,n1) b(n1,1) complex
expression Objective
Objective(k) = square_pos(norm((K*x(:,1)-b(:,1)),2)) + lambda^2 * square_pos( norm(x(:,1),2)); 
minimize(sum(Objective))   
cvx_end
toc

figure();
subplot(2,1,1);
plot(abs(Ex))
title('(x');
subplot(2,1,2); 
plot(E)
title('(E)');

(shine_wade) #9

maybe the parameters are a bit difference ,but the error make me confused


(Mark L. Stone) #10

The last of these lines should be first

m1=size(K,1);  
n1=size(K,2);  
K=Tw.*(aep-az)*ep0*epr;

The statement
variables x(n1,1) K(m1,n1) b(n1,1) complex
is not allowed because CVX thinks you are declaring
complex
as a CVX variable, but that is not allowed because complex is a reserved keyword. You need to use
variable
rather than
variables
if you wish to declare them as complex.
http://cvxr.com/cvx/doc/basics.html#variables

The one limitation of the variables command is that it cannot declare complex, integer, or structured variables. These must be declared one at a time, using the singular variable command.

Using
Objective(k()
is wrong, but easily fixed.

But your big difficulty is having
K*x(:,1)
because K and x are both declared as CVX variables. That violates CVX’s rules and by itself is non-convex…However, i believe you want to use K as previously defined, and not declare it in CVX, which would overrides its previous value. Then your program will run. It’s still running on my machine, but I’ll post this now.

cvx_begin
variable x(n1,1) complex
variable b(n1,1) complex
Objective = 0;
for joe_blow = 1:n1
  Objective = Objective + square_pos(norm((K*x(joe_blow,1)-b(joe_blow,1)),2)) + lambda^2 * square_pos( norm(x(joe_blow,1),2)); 
end
minimize(Objective)   
cvx_end

Is your problem also supposed to have (linear equality) constraints?


(shine_wade) #11

thankyou , Mark ,very appieciate ,you have solved my first error:Disciplined convex programming error:
Invalid quadratic form: not a scalar. but after I run the changed code for a long time,it show :*In_Incorrect use cvx/subsref (line 19)The index is beyond the matrix dimension..To speed up, I use vectors instead of loops, but the error also come:Incorrect use cvx/subsref (line 19)The index is beyond the matrix dimension.**

the code show as below:

tic;
close all
format short
thmcd=0.22; 	
dens=910;   	
c=1900;     		
A=1.54*10^-4;  	    
yita=0.1;   		
q=0.005;    		
tcth=0.2;   		
D=1.2724*1e-7;
e1=10*1e6;
d=9.8*1e-6;      	
multi=1;
D=D*multi;   				
StartPosition=9.3e-8*exp(0.8e5*d)+1.2e-6;

tc=(d)^2/(pi^2*D);  			
T0=yita*q/(c*dens*A*d);   
az=1.35*10^-4;  		
aep=2.45*az;   		
ep0=8.8542*10^-12;  	
epr=2.2;    			
ts=1e-7;       		
tstart=1e-7;   		
tstop=0.005;   		
t=tstart:ts:tstop;
z=linspace(0,d,100); 

fmin=1/tstop;
fmax=4*10^5;    
f=logspace(log10(fmin),log10(fmax),400);
E1=linspace(0,10,25);
E2=linspace(10,10,length(z)-56);
E3=linspace(10,0,31);
E=[E1 E2 E3];
G=(aep-az)*ep0*epr*E;
n=1:100;
T=T0*ones(length(z),1)*exp(-t/tcth)...             	 	
    +2*T0*cos(pi*z'*n/d)*exp(-(n.^2)'*t/tc);
dTdt=-T0/tcth*ones(length(z),1)*exp(-t/tcth)...     		
    +2*T0*(ones(length(z),1)*(-n.^2/tc)).*cos(pi*z'*n/d)*exp(-(n.^2)'*t/tc);
Gtemp=zeros(1,length(z)-1);
Tt1=zeros(length(z)-1,1);
dTdttemp=zeros(length(z)-1,length(t));
for k=1:length(z)-1
    Gtemp(k)=(G(k)+G(k+1))/2;
    Tt1(k)=(T(k,1)+T(k+1,1))/2;
    dTdttemp(k,:)=(dTdt(k,:)+dTdt(k+1,:))/2;
end

If=A/(length(z)-1)*Gtemp*dTdttemp;
Qt1f=A/(length(z)-1)*Gtemp*Tt1;
nf=length(f);
Yf=ts*If*exp(-1i*2*pi*t'*f)+Qt1f; 
Yfr=real(Yf);
Yfi=imag(Yf);
w=2*pi*f;
k=(1+1i)*(sqrt(w./(2*D)))';
ks=k./sinh(k.*d);
Tw0=yita*q/(c*dens*A);  
Tw=Tw0*zeros(size(f,2),size(z,2))+ks.*cosh(k.*(d-z));

%Kx=b
m1=size(K,1);  
n1=size(K,2);  
K=Tw.*(aep-az)*ep0*epr;
b=Yf';
%K;   %m1 x n1 matrix    
lambda = 1.045929285995968e-09;

%Use cvx
tic
cvx_begin
variable x(n1,1) complex
variable b(m1,1) complex
joe_blow = 1:m1;
Objective = zeros(n1,1)+square_pos(norm((K*x(joe_blow,1)-b(:,1)),2)) + lambda^2 * square_pos( norm(x(joe_blow,1),2)); 
% Objective = 0;
% for joe_blow = 1:m1
%   Objective = Objective + square_pos(norm((K*x(joe_blow,1)-b(joe_blow,1)),2)) + lambda^2 * square_pos( norm(x(joe_blow,1),2)); 
% end
minimize(Objective)  
cvx_end
toc

figure();
subplot(2,1,1);
plot(abs(Ex))
title('(x');
subplot(2,1,2); 
plot(E)
title('(E)');

(Mark L. Stone) #12

You should have
joe_blow = 1:n1;
not
joe_blow = 1:m1;

Also, you need to sum Objective.

Here is the corrected program and output. I leave it to you to determine whether this correctly implements the model you wish to solve.

tic
cvx_begin
variable x(n1,1) complex
variable b(m1,1) complex
joe_blow = 1:n1;
Objective = zeros(n1,1)+square_pos(norm((K*x(joe_blow,1)-b(:,1)),2)) + lambda^2 * square_pos( norm(x(joe_blow,1),2)); 
% Objective = 0;
% for joe_blow = 1:m1
%   Objective = Objective + square_pos(norm((K*x(joe_blow,1)-b(joe_blow,1)),2)) + lambda^2 * square_pos( norm(x(joe_blow,1),2)); 
% end
minimize(sum(Objective))
cvx_end
toc
Calling SDPT3 4.0: 1012 variables, 1008 equality constraints
   For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------

 num. of constraints = 1008
 dim. of sdp    var  =  4,   num. of sdp  blk  =  2
 dim. of socp   var  = 1002,   num. of socp blk  =  2
 dim. of linear var  =  4
*******************************************************************
   SDPT3: Infeasible path-following algorithms
*******************************************************************
 version  predcorr  gam  expon  scale_data
   HKM      1      0.000   1        0    
it pstep dstep pinfeas dinfeas  gap      prim-obj      dual-obj    cputime
-------------------------------------------------------------------
 0|0.000|0.000|5.2e+01|1.7e+01|1.2e+05| 3.030000e+02  0.000000e+00| 0:0:00| spchol  1  1 
 1|0.631|0.613|1.9e+01|6.9e+00|4.5e+04| 4.550192e+02 -1.372928e+03| 0:0:00| spchol  1  1 
 2|0.703|0.731|5.7e+00|1.9e+00|1.4e+04| 1.193615e+03 -2.747421e+03| 0:0:00| spchol  1  1 
 3|0.793|0.874|1.2e+00|2.6e-01|2.8e+03| 9.908874e+02 -6.322516e+02| 0:0:00| spchol  1  1 
 4|0.652|0.615|4.1e-01|1.1e-01|9.8e+02| 3.874594e+02 -2.783146e+02| 0:0:00| spchol  1  1 
 5|0.482|0.752|2.1e-01|2.8e-02|4.1e+02| 1.867971e+02 -1.134500e+02| 0:0:01| spchol  1  1 
 6|0.498|0.930|1.1e-01|2.2e-03|1.6e+02| 8.046219e+01 -3.996481e+01| 0:0:01| spchol  1  1 
 7|0.552|1.000|4.8e-02|2.4e-05|6.0e+01| 2.923157e+01 -1.324480e+01| 0:0:01| spchol  1  1 
 8|0.596|0.875|1.9e-02|9.6e-03|2.3e+01| 1.009677e+01 -4.411257e+00| 0:0:01| spchol  1  1 
 9|0.769|0.662|4.5e-03|7.1e-03|6.7e+00| 2.835247e+00 -1.237396e+00| 0:0:01| spchol  1  1 
10|0.779|0.660|9.9e-04|3.3e-03|1.9e+00| 8.409754e-01 -3.051821e-01| 0:0:01| spchol  1  1 
11|0.774|0.637|2.2e-04|1.4e-03|5.4e-01| 2.548194e-01 -7.518300e-02| 0:0:01| spchol  1  1 
12|0.756|0.623|5.5e-05|5.7e-04|1.7e-01| 7.935999e-02 -1.873660e-02| 0:0:01| spchol  1  1 
13|0.732|0.618|1.5e-05|2.3e-04|5.5e-02| 2.581796e-02 -4.453933e-03| 0:0:01| spchol  1  1 
14|0.725|0.612|4.0e-06|9.2e-05|1.8e-02| 8.605029e-03 -7.554578e-04| 0:0:01| spchol  1  1 
15|0.732|0.605|1.1e-06|3.7e-05|6.1e-03| 2.878326e-03  9.405693e-05| 0:0:01| spchol  1  1 
16|0.748|0.597|2.7e-07|1.5e-05|2.0e-03| 9.551918e-04  1.981771e-04| 0:0:01| spchol  1  1 
17|0.775|0.587|6.1e-08|6.3e-06|6.7e-04| 3.110299e-04  1.479872e-04| 0:0:01| spchol  1  1 
18|0.831|0.572|1.0e-08|2.7e-06|2.2e-04| 9.678154e-05  9.181818e-05| 0:0:01| spchol  1  1 
19|0.911|0.547|9.2e-10|1.2e-06|7.2e-05| 3.182792e-05  5.494112e-05| 0:0:01| spchol  1  1 
20|0.834|0.605|1.5e-10|4.9e-07|2.2e-05| 9.692306e-06  2.471075e-05| 0:0:01| spchol  1  1 
21|0.887|0.547|1.7e-11|2.2e-07|7.8e-06| 3.715018e-06  1.276063e-05| 0:0:01| spchol  1  1 
22|0.465|0.761|9.3e-12|5.3e-08|3.1e-06| 1.823485e-06  2.760187e-06| 0:0:01| spchol  1  1 
23|0.497|1.000|4.7e-12|1.9e-12|1.3e-06| 5.676709e-07 -7.151980e-07| 0:0:02| spchol  1  1 
24|0.488|0.724|2.4e-12|1.5e-12|4.8e-07| 3.157197e-07 -1.631266e-07| 0:0:02| spchol  1  1 
25|0.483|1.000|1.2e-12|1.0e-12|2.1e-07| 1.210256e-07 -8.927485e-08| 0:0:02| spchol  1  1 
26|0.480|0.928|6.4e-13|1.1e-12|8.0e-08| 5.404697e-08 -2.411331e-08| 0:0:02| spchol  1  1 
27|0.478|1.000|3.4e-13|1.0e-12|3.8e-08| 2.028837e-08 -1.684388e-08| 0:0:02| spchol  1  1 
28|0.476|1.000|1.8e-13|1.0e-12|1.4e-08| 7.954901e-09 -5.717309e-09| 0:0:02|
  stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
 number of iterations   = 28
 primal objective value =  7.95490079e-09
 dual   objective value = -5.71730884e-09
 gap := trace(XZ)       = 1.42e-08
 relative gap           = 1.42e-08
 actual relative gap    = 1.37e-08
 rel. primal infeas (scaled problem)   = 1.76e-13
 rel. dual     "        "       "      = 1.00e-12
 rel. primal infeas (unscaled problem) = 0.00e+00
 rel. dual     "        "       "      = 0.00e+00
 norm(X), norm(y), norm(Z) = 1.0e+02, 6.4e+02, 6.4e+02
 norm(A), norm(b), norm(C) = 3.3e+01, 1.0e+02, 2.4e+00
 Total CPU time (secs)  = 1.80  
 CPU time per iteration = 0.06  
 termination code       =  0
 DIMACS: 1.8e-13  0.0e+00  1.2e-12  0.0e+00  1.4e-08  1.4e-08
-------------------------------------------------------------------
 
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +5.71731e-09
 
Elapsed time is 2.482381 seconds.