# Ideas for BP
versioninfo()
# 生成数据
using LinearAlgebra
using Distributions
using ForwardDiff
# 导入随机数生成包
using Random
# 系数矩阵的行数和列数
m = 100
n = 1024
# 信号的稀疏度
s = 10
# 稀疏信号的构造
xo = zeros(n,1);
pn = randperm(n);
inds = pn[1:s];
xo[inds] = randn(s,1);
# 系数矩阵的构造
A = randn(m,n);
ε = 1e-6
A .-= ε*mean(A, dims = 1); # this apparently is good practice, better conditioning
A ./= sqrt.(sum(abs2, A, dims = 1));
# 生成观察数据
b = A*xo;
println("数据构造完毕,系数矩阵A和观察数据b\n")
# 给出原始数据的一些结果
println("L1 of xo : \n",norm(xo,1))
# 绘制原始数据的图像
using Plots
pxo = plot(1:n,xo,linewidth = 2,label="Original",shape=:circle,color=:red)
# savefig(pxo, "BPorignal.pdf")
# 直接求解非线性规划问题,即BP
using JuMP
# 导入求解器
using Ipopt
# 构造一个模型对象
BP = Model(Ipopt.Optimizer)
set_optimizer_attribute(BP, "print_level", 4)
# 设置变量
@variable(BP,bpx[1:n]);
# 设置目标函数,直接求解L1范数
# @NLobjective(BP,Min,sum(abs(bpx[i]) for i =1:n))
# 求解L1/L2 或者是 L1/Linf
l1oL2(x::Vector) = norm(x,1)/(norm(x,2)+0.000001)
l1olinf(x::Vector) = norm(x,1)/(norm(x,Inf)+0.0000001)
# Register
register(BP, :wjL1L2, 1, l1oL2; autodiff = true)
register(BP, :wjL1Linf, 1, l1olinf; autodiff = true)
# 构造目标函数
@NLobjective(BP,Min,wjL1L2(bpx))
# 设置约束条件
@constraint(BP, b .== A*bpx);
# 打印模型
# print(BP)
# 求解
resBP = JuMP.optimize!(BP)
solution_summary(BP)
# 输出解
BPx = JuMP.value.(bpx)
# 绘制求解的解的图示
BPx = JuMP.value.(bpx)
# 绘制求解的解的图示
psolBP = plot(1:n,BPx,linewidth = 2,label="BP_Ipopt",shape=:star5,color=:green)
问题
目标函数我想定义为l1范数和l2范数之比,但是就是不能通过!!!!!!!!!
现在忽然发现,Julia也没有想想中的那么好,之前看到一个博客( 为什么我不再推荐你用Julia? (qq.com)),或许是有原因的.话又说回来了,真的有那么好的,我不信不会有很多人放弃一些语言转投Julia,但并没有出现这种,所以,她或许没有那么好!而我也隐约有个不好的预感:突然有一天,julia被放弃啦!