Skip to content

Commit

Permalink
More efficient and less intrusive quadform (#444)
Browse files Browse the repository at this point in the history
* small change to make it work

* Update src/atoms/second_order_cone/quadform.jl

Co-authored-by: Eric Hanson <5846501+ericphanson@users.noreply.github.com>

* Update src/atoms/second_order_cone/quadform.jl

Co-authored-by: Eric Hanson <5846501+ericphanson@users.noreply.github.com>

* Adding a test

* Update src/problem_depot/problems/socp.jl

Co-authored-by: Eric Hanson <5846501+ericphanson@users.noreply.github.com>

Co-authored-by: Eric Hanson <5846501+ericphanson@users.noreply.github.com>
  • Loading branch information
lrnv and ericphanson committed Jul 5, 2021
1 parent feacb66 commit 73aec9c
Show file tree
Hide file tree
Showing 2 changed files with 18 additions and 7 deletions.
12 changes: 5 additions & 7 deletions src/atoms/second_order_cone/quadform.jl
Original file line number Diff line number Diff line change
Expand Up @@ -6,20 +6,18 @@ function quadform(x::AbstractExpr, A::Value)
if length(size(A)) != 2 || size(A, 1) != size(A, 2)
error("Quadratic form only takes square matrices")
end
if !issymmetric(A)
error("Quadratic form only defined for symmetric matrices")
if !ishermitian(A)
error("Quadratic form only defined for Hermitian matrices")
end
V = eigvals(Symmetric(Matrix(A)))

if all(V .>= 0)
if isposdef(A)
factor = 1
elseif all(V .<= 0)
elseif isposdef(.-A)
factor = -1
else
error("Quadratic forms supported only for semidefinite matrices")
end

P = real(sqrt(Matrix(factor * A)))
P = sqrt(Hermitian(factor * A))
return factor * square(norm2(P * x))
end

Expand Down
13 changes: 13 additions & 0 deletions src/problem_depot/problems/socp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -267,6 +267,19 @@ end
@test evaluate(H) Hval atol=atol rtol=rtol
end

# https://github.com/jump-dev/Convex.jl/pull/444
x = Variable(3)
a,b,c,d,e,f = rand(6)
M = [2 a-b*im c-d*im;
a+b*im 2 e-f*im
c+d*im e+f*im 2]
y = rand(3)
p = minimize(quadform(x-y,M); numeric_type = T)
handle_problem!(p)
if test
@test evaluate(x) y atol=atol rtol=rtol
end

end

@add_problem socp function socp_huber_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
Expand Down

0 comments on commit 73aec9c

Please sign in to comment.