diff --git a/Project.toml b/Project.toml
index b5269b5e..a12b1a1e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -20,13 +20,13 @@ Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 Sobol = "ed01d8cd-4d21-5b2a-85b4-cc3bdc58bad4"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
-Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [compat]
 ADNLPModels = "0.1, 0.2, 0.3"
@@ -55,6 +55,8 @@ Zygote = "0.6"
 julia = "1"
 
 [extras]
+ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 Hyperopt = "93e5fe13-2215-51db-baaf-2e9a34fb2712"
 Juniper = "2ddba703-00a4-53a7-87a5-e8b9971dde84"
@@ -65,4 +67,4 @@ SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["FiniteDifferences", "Hyperopt", "Juniper", "NLopt", "Pavito", "Percival", "SafeTestsets", "Test"]
+test = ["FiniteDifferences", "Hyperopt", "Juniper", "NLopt", "Pavito", "Percival", "SafeTestsets", "Distributions", "Test", "ChainRulesTestUtils"]
diff --git a/src/Nonconvex.jl b/src/Nonconvex.jl
index a32db59d..6adc5d96 100644
--- a/src/Nonconvex.jl
+++ b/src/Nonconvex.jl
@@ -9,6 +9,8 @@ export  Model,
         set_objective!,
         add_ineq_constraint!,
         add_eq_constraint!,
+        add_sd_constraint!,
+        decompress_symmetric,
         getmin,
         getmax,
         isinteger,
@@ -32,6 +34,7 @@ export  Model,
         PavitoIpoptCbcAlg,
         HyperoptAlg,
         BayesOptAlg,
+        SDPBarrierAlg,
         MTSAlg,
         LS1Alg,
         KKTCriteria,
@@ -47,6 +50,7 @@ export  Model,
         PavitoIpoptCbcOptions,
         HyperoptOptions,
         BayesOptOptions,
+        SDPBarrierOptions,
         MTSOptions,
         LS1Options,
         Tolerance,
@@ -126,6 +130,9 @@ include("algorithms/mts.jl")
 
 include("algorithms/deflation.jl")
 
+# Semi-definite programming
+include("algorithms/sdp.jl")
+
 # Wrappers
 
 include("wrappers/moi.jl")
diff --git a/src/algorithms/common.jl b/src/algorithms/common.jl
index 2ef51bb6..376b3482 100644
--- a/src/algorithms/common.jl
+++ b/src/algorithms/common.jl
@@ -212,6 +212,7 @@ end
 
 abstract type AbstractModel end
 
+
 """
 ```
 optimize(
@@ -224,6 +225,7 @@ optimize(
     callback::Function = plot_trace ? LazyPlottingCallback() : NoCallback(),
     kwargs...,
 )
+
 ```
 
 Optimizes `model` using the algorithm `optimizer`, e.g. an instance of [`MMA87`](@ref) or [`MMA02`](@ref). `x0` is the initial solution. The keyword arguments are:
@@ -234,7 +236,14 @@ Optimizes `model` using the algorithm `optimizer`, e.g. an instance of [`MMA87`]
 
  The details of the MMA optimization algorithms can be found in the original [1987 MMA paper](https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620240207) and the [2002 paper](https://epubs.siam.org/doi/abs/10.1137/S1052623499362822).
 """
+function _optimize_precheck(model::AbstractModel, optimizer::AbstractOptimizer, args...; options=nothing, kwargs...)
+    if (length(model.sd_constraints.fs) != 0) && !(options isa SDPBarrierOptions)
+        @warn "Input `sd_constraints` will be ignored due to not using `SDPBarrierOptions`. "
+    end
+end
+
 function optimize(model::AbstractModel, optimizer::AbstractOptimizer, x0, args...; kwargs...)
+    _optimize_precheck(model, optimizer, x0, args...; kwargs...)
     _model, _x0, unflatten = tovecmodel(model, x0)
     r = optimize(_model, optimizer, _x0, args...; kwargs...)
     return @set r.minimizer = unflatten(r.minimizer)
@@ -244,11 +253,13 @@ end
  optimize without x0
 """
  function optimize(model::AbstractModel, optimizer::AbstractOptimizer, args...; kwargs...)
+    _optimize_precheck(model, optimizer, args...; kwargs...)
     _model, _, unflatten = tovecmodel(model)
     r = optimize(_model, optimizer, args...; kwargs...)
     return @set r.minimizer = unflatten(r.minimizer)
 end
 
+
 """
 Clamp the box constraint and evaluate objective function at point x
 """
diff --git a/src/algorithms/mma_algorithm.jl b/src/algorithms/mma_algorithm.jl
index de22c0aa..18935fcd 100755
--- a/src/algorithms/mma_algorithm.jl
+++ b/src/algorithms/mma_algorithm.jl
@@ -235,6 +235,7 @@ function optimize!(workspace::MMAWorkspace)
 
             # Evaluates the exact objective and constraints and their gradients at the optimal x
             optimalx = getoptimalx(dualmodel)
+
             fg, ∇fg = value_jacobian(
                 approxmodel.objective_ineq_constraints,
                 optimalx,
@@ -245,7 +246,6 @@ function optimize!(workspace::MMAWorkspace)
             # Reduces the chance of overflow or underflow in λs 
             # Helps with the convergence
             auto_scale && iter < 5 && scaleobjective!(model, fg, ∇fg, approxfg, λ)
-            #iter < 10 && scaleobjective!(model, fg, ∇fg, approxfg, λ)
 
             # Check if the approximation is an upper bound at the current x.
             # If not, increase ρ. Only used in MMA02
diff --git a/src/algorithms/sdp.jl b/src/algorithms/sdp.jl
new file mode 100644
index 00000000..d36c7e17
--- /dev/null
+++ b/src/algorithms/sdp.jl
@@ -0,0 +1,137 @@
+# Semidefinite programming
+
+# Decompress
+function lowertriangind(mat::Matrix)
+    indices = [i for i in CartesianIndices(mat) if i[1]>i[2]]
+    return LinearIndices(mat)[indices]
+end
+
+function rearrange_x(x_L::AbstractVector, x_D::AbstractVector)
+    mat_dim = length(x_D)
+    L = zeros(mat_dim, mat_dim)
+    L[lowertriangind(L)] .= x_L
+    D = diagm(x_D)
+    return (L, D)
+end
+
+function decompress_symmetric(L::Matrix, D::Matrix)
+    L + D + L'
+end
+
+function decompress_symmetric(x_L::AbstractArray, x_D::AbstractArray)
+    L, D = rearrange_x(x_L, x_D)
+    return decompress_symmetric(L, D)
+end
+
+function ChainRulesCore.rrule(::typeof(rearrange_x), x_L::AbstractVector, x_D::AbstractVector)
+    function pullback((ΔL, ΔD))
+        Δx_L = ΔL[lowertriangind(ΔL)]
+        Δx_D = diag(ΔD)
+        NoTangent(), Δx_L, Δx_D 
+    end
+    return rearrange_x(x_L, x_D), pullback
+end
+
+
+# Optimizer (algorithm)
+struct SDPBarrierAlg{Alg <: AbstractOptimizer} <: AbstractOptimizer
+    sub_alg::Alg
+end
+function SDPBarrierAlg(;sub_alg)
+    return SDPBarrierAlg(sub_alg)
+end
+
+# Options
+@params mutable struct SDPBarrierOptions
+    # Dimension of objective matrix
+    # Hyperparameters 
+    # Initial value of `c` in barrier method: 
+    # `Real` for using same value for all `sd_constraints`, `AbstractArray` for assign them respectively
+    c_init::Union{Real, AbstractArray}
+    # Decrease rate of `c` for every epoch, same as above
+    c_decr::Union{Real, AbstractArray}
+    n_iter::Int
+    # sub_option to solve (in)equality constraints
+    sub_options
+    # Keep all results or not
+    keep_all::Bool
+end
+function SDPBarrierOptions(c_init, c_decr, n_iter; sub_options, keep_all=true)
+    @assert 0 < c_decr < 1 "c_decr should be between 0 and 1. "
+    @assert c_init > 0 "c_init shoule be larger than 0. "
+    SDPBarrierOptions(c_init, c_decr, n_iter, sub_options, keep_all)
+end
+function SDPBarrierOptions(;sub_options, c_init=1.0, c_decr=0.1, n_iter=10, keep_all=true)
+    SDPBarrierOptions(c_init, c_decr, n_iter, sub_options=sub_options, keep_all=keep_all)
+end
+
+# Result
+@params struct SDPResult
+    minimum
+    minimizer
+    results
+    optimal_ind
+end
+
+@params struct SDPWorkspace <: Workspace
+    model::VecModel
+    x0::AbstractVector
+    options::SDPBarrierOptions
+    sub_alg::AbstractOptimizer
+end
+
+function Workspace(model::VecModel, optimizer::SDPBarrierAlg, x0, args...; options, kwargs...,)
+    @unpack c_init, c_decr = options
+    for c in model.sd_constraints.fs
+        @assert isposdef(c(x0)) "Initial matrix should be positive definite. "
+    end
+    if c_init isa AbstractArray
+        @assert length(model.sd_constraints.fs) == length(c_init) "c_init should be same length with number of `sd_constraints` when using array. "
+    end
+    if c_decr isa AbstractArray
+        @assert length(model.sd_constraints.fs) == length(c_decr) "c_decr should be same length with number of `sd_constraints` when using array. "
+    end
+    if c_init isa AbstractArray && c_decr isa AbstractArray
+        @assert length(c_init) == length(c_decr) "c_decr should be same length with c_init. "
+    end
+    return SDPWorkspace(model, copy(x0), options, optimizer.sub_alg)
+end
+
+function sd_objective(objective0, sd_constraints, c::AbstractArray)
+    function _objective(args)
+        target = objective0(args)
+        barrier = sum(c .* -logdet.(map(f -> f(args), sd_constraints.fs)))
+        return target + barrier
+    end
+    return _objective
+end
+
+function to_barrier(model, c::AbstractArray)
+    sd_constraints, objective0 = model.sd_constraints, model.objective
+    _model = set_objective(model, sd_objective(objective0, sd_constraints, c))
+    return _model
+end
+
+function optimize!(workspace::SDPWorkspace)
+    @unpack model, x0, options, sub_alg = workspace
+    @unpack c_init, c_decr, n_iter, sub_options, keep_all = options
+    objective0 = model.objective
+    x = copy(x0)
+    c = c_init isa Real ? ([c_init for _ in 1:length(model.sd_constraints.fs)]) : c_init
+    results = []
+    for _ in 1:n_iter
+        model_i = to_barrier(model, c)
+        result_i = optimize(model_i, sub_alg, x, options = sub_options)
+        minimizer_i = result_i.minimizer
+        push!(results, (objective0(minimizer_i), minimizer_i))
+        c = c .* c_decr
+        x = copy(minimizer_i)
+    end
+    optimal_ind = argmin(first.(results))
+    minimum, minimizer = results[optimal_ind]
+    if keep_all
+        return SDPResult(minimum, minimizer, results, optimal_ind)
+    else
+        return SDPResult(minimum, minimizer, nothing, nothing)
+    end
+end
diff --git a/src/functions/functions.jl b/src/functions/functions.jl
index 078609ba..78f48861 100755
--- a/src/functions/functions.jl
+++ b/src/functions/functions.jl
@@ -273,3 +273,25 @@ getdim(c::EqConstraint) = c.dim
 Returns the function wrapped in `f`.
 """
 getfunction(f::EqConstraint) = f.f
+
+
+"""
+
+Used in semidefinite programming
+"""
+@params struct SDConstraint <: AbstractConstraint
+    f::Function
+    dim::Int
+end
+
+function (c::SDConstraint)(args...; kwargs...)
+    out = c.f(args...; kwargs...)
+    # println(length(out))
+    # println(getdim(c))
+    @assert length(out) == getdim(c)^2
+    return out
+end
+
+getdim(c::SDConstraint) = c.dim
+
+getfunction(c::SDConstraint) = c.f
diff --git a/src/models/dict_model.jl b/src/models/dict_model.jl
index 9f4cf246..a74f97ad 100644
--- a/src/models/dict_model.jl
+++ b/src/models/dict_model.jl
@@ -2,6 +2,7 @@ mutable struct DictModel <: AbstractModel
     objective::Union{Nothing, Objective}
     eq_constraints::VectorOfFunctions
     ineq_constraints::VectorOfFunctions
+    sd_constraints::VectorOfFunctions
     box_min::OrderedDict
     box_max::OrderedDict
     init::OrderedDict
@@ -12,6 +13,7 @@ function DictModel(f = nothing)
     return DictModel(
         Objective(f), VectorOfFunctions(EqConstraint[]),
         VectorOfFunctions(IneqConstraint[]),
+        VectorOfFunctions(SDConstraint[]),
         OrderedDict(), OrderedDict(),
         OrderedDict(), OrderedDict()
     )
diff --git a/src/models/model.jl b/src/models/model.jl
index aac782c9..238b3163 100755
--- a/src/models/model.jl
+++ b/src/models/model.jl
@@ -15,11 +15,13 @@ The `Model` structs stores information about the nonlinear optimization problem.
 - `objective`: the objective function of the problem of type [`Objective`](@ref).
 - `eq_constraints`: the equality constraints of the problem of type [`VectorOfFunctions`](@ref). Each function in `ineq_constraints` can be an instance of [`IneqConstraint`](@ref) or [`AbstractFunction`](@ref). If the function is not an `IneqConstraint`, its right-hand-side bound is assumed to be 0.
 - `ineq_constraints`: the inequality constraints of the problem of type [`VectorOfFunctions`](@ref). Each function in `ineq_constraints` can be an instance of [`IneqConstraint`](@ref) or [`AbstractFunction`](@ref). If the function is not an `IneqConstraint`, its right-hand-side bound is assumed to be 0.
+- `sd_constraints`: a list of matrix-valued functions which must be positive semidefinite at any feasible solution.
 """
 mutable struct Model{Tv <: AbstractVector} <: AbstractModel
     objective::Union{Nothing, Objective}
     eq_constraints::VectorOfFunctions
     ineq_constraints::VectorOfFunctions
+    sd_constraints::VectorOfFunctions
     box_min::Tv
     box_max::Tv
     init::Tv
@@ -41,7 +43,12 @@ Constructs an empty model with objective function `f` and decision variable valu
 """
 Model(::Type{T}, f::Function) where {T} = Model(T, Objective(f))
 function Model(::Type{T}, obj::Union{Nothing, Objective}) where {T}
-    return Model(obj, VectorOfFunctions(EqConstraint[]), VectorOfFunctions(IneqConstraint[]), T[], T[], T[], falses(0))
+    return Model(obj, 
+    VectorOfFunctions(EqConstraint[]), 
+    VectorOfFunctions(IneqConstraint[]), 
+    VectorOfFunctions(SDConstraint[]), 
+    T[], T[], T[], 
+    falses(0))
 end
 
 getobjective(m::AbstractModel) = m.objective
@@ -139,6 +146,11 @@ end
 
 Sets the objective of the moodel `m` to the function `f`. `f` must return a scalar.
 """
+function set_objective(m::AbstractModel, f::Function; kwargs...)
+    @set m.objective = Objective(f; kwargs...)
+    return m
+end
+
 function set_objective!(m::AbstractModel, f::Function; kwargs...)
     m.objective = Objective(f; kwargs...)
     return m
@@ -173,3 +185,11 @@ function add_eq_constraint!(m::AbstractModel, fs::Vector{<:EqConstraint})
     append!(m.eq_constraints.fs, fs)
     return m
 end
+
+function add_sd_constraint!(m::AbstractModel, f::Function, dim=size(f(getinit(m)), 1))
+    return add_sd_constraint!(m, SDConstraint(f, dim))
+end
+function add_sd_constraint!(m::AbstractModel, sd_constraint::AbstractFunction)
+    push!(m.sd_constraints.fs, sd_constraint)
+    return m
+end
diff --git a/src/models/vec_model.jl b/src/models/vec_model.jl
index 0008c8fe..f11522ed 100644
--- a/src/models/vec_model.jl
+++ b/src/models/vec_model.jl
@@ -2,13 +2,21 @@ mutable struct VecModel{Tv <: AbstractVector} <: AbstractModel
     objective::Union{Nothing, Objective}
     eq_constraints::VectorOfFunctions
     ineq_constraints::VectorOfFunctions
+    sd_constraints::VectorOfFunctions
     box_min::Tv
     box_max::Tv
     init::Tv
     integer::BitVector
 end
-function VecModel(objective, eq_constraints, ineq_constraints, box_min, box_max, init, integer)
-    VecModel(objective, eq_constraints, ineq_constraints, box_min, box_max, init, integer, length(box_min))
+function VecModel(
+        objective::Union{Nothing, Objective}, 
+        eq_constraints::VectorOfFunctions, 
+        ineq_constraints::VectorOfFunctions, 
+        box_min, 
+        box_max, 
+        init, 
+        integer::BitVector)
+    return VecModel(objective, eq_constraints, ineq_constraints, VectorOfFunctions(SDConstraint[]), box_min, box_max, init, integer)
 end
 
 function isfeasible(model::VecModel, x::AbstractVector; ctol = 1e-4)
@@ -57,6 +65,9 @@ function set_objective_multiple!(model::VecModel, m)
     return model
 end
 
+"""
+Generic `optimize` for VecModel
+"""
 function optimize(model::VecModel, optimizer::AbstractOptimizer, x0, args...; kwargs...)
     workspace = Workspace(model, optimizer, copy(x0), args...; kwargs...)
     return optimize!(workspace)
@@ -74,16 +85,27 @@ function tovecmodel(m::AbstractModel, x0 = getmin(m))
     v, _unflatten = flatten(x0)
     unflatten = Unflatten(x0, _unflatten)
     return VecModel(
+        # objective
         Objective(x -> m.objective(unflatten(x)), flags = m.objective.flags),
+        # eq_constraints
         length(m.eq_constraints.fs) != 0 ? VectorOfFunctions(map(m.eq_constraints.fs) do c
             EqConstraint(x -> maybeflatten(c.f(unflatten(x)))[1], maybeflatten(c.rhs)[1], c.dim, c.flags)
         end) : VectorOfFunctions(EqConstraint[]),
+        # ineq_constraints
         length(m.ineq_constraints.fs) != 0 ? VectorOfFunctions(map(m.ineq_constraints.fs) do c
             IneqConstraint(x -> maybeflatten(c.f(unflatten(x)))[1], maybeflatten(c.rhs)[1], c.dim, c.flags)
         end) : VectorOfFunctions(IneqConstraint[]),
+        # sd_constraints
+        length(m.sd_constraints.fs) != 0 ? VectorOfFunctions(map(m.sd_constraints.fs) do c
+            SDConstraint(x -> c.f(unflatten(x)), c.dim)
+        end) : VectorOfFunctions(SDConstraint[]),
+        # box_min
         float.(flatten(m.box_min)[1]),
+        # box_max
         float.(flatten(m.box_max)[1]),
+        # init
         float.(flatten(m.init)[1]),
+        # integer
         convert(BitVector, flatten(m.integer)[1]),
     ), float.(v), unflatten
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 53006fa2..7dfe3f80 100755
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -17,3 +17,4 @@ end
 @safetestset "Bayesian optimization" begin include("bayesian.jl") end
 @safetestset "Multiple trajectory search" begin include("mts.jl") end
 @safetestset "Deflation" begin include("deflation.jl") end
+@safetestset "Semidefinite programming" begin include("sdp.jl") end
diff --git a/test/sdp.jl b/test/sdp.jl
new file mode 100644
index 00000000..c30df71f
--- /dev/null
+++ b/test/sdp.jl
@@ -0,0 +1,146 @@
+# Test cases of semidefinite programming
+using Nonconvex, LinearAlgebra, Test, Zygote, FiniteDifferences, Distributions
+using ChainRulesTestUtils, Random
+Random.seed!(1)
+
+const FDM = FiniteDifferences
+
+# Test setting
+mat_dim = 2
+mat_length = mat_dim^2
+n_sample = 300
+
+@testset "Test matrix to vectors function" begin
+    test_rrule(Nonconvex.rearrange_x, rand((mat_dim^2-mat_dim)÷2), rand(mat_dim))
+end
+
+function random_psd_mat(mat_dim)
+    _mat = randn(mat_dim, mat_dim)
+    return _mat' * _mat
+end
+function sd_constraint((x_L, x_D))
+    decompress_symmetric(x_L, x_D)
+end
+
+@testset "Single semi-definite constraint" begin
+    # Randomize groundtruth
+    μ, Σ = randn(mat_dim), random_psd_mat(mat_dim)
+
+    # Generate samples
+    samples = rand(MvNormal(μ, Σ), n_sample)
+
+    # Objective function
+    function f((x_L, x_D))
+        -loglikelihood(MvNormal(μ, decompress_symmetric(x_L, x_D)), samples)
+    end
+
+    model = Model(f)
+
+    mat_x0 = random_psd_mat(mat_dim)
+    x0 = [mat_x0[Nonconvex.lowertriangind(mat_x0)], diag(mat_x0)]
+    lbs = [fill(-Inf, length(x0[1])), zeros(length(x0[2]))]
+    ubs = [fill(Inf, length(x0[1])), fill(Inf, length(x0[2]))]
+    addvar!(model, lbs, ubs)
+
+    add_sd_constraint!(model, sd_constraint)
+
+    alg = SDPBarrierAlg(sub_alg=IpoptAlg())
+    options = SDPBarrierOptions(sub_options=IpoptOptions(max_iter=200))
+    result = optimize(model, alg, x0, options = options)
+
+    minimum, minimizer, optimal_ind = result.minimum, result.minimizer, result.optimal_ind
+    _Σ = sd_constraint(minimizer)
+
+    println("result: \n $result")
+
+    println("minimum: \n $minimum")
+    println("minimizer: \n $minimizer")
+    println("_Σ: \n $(_Σ)")
+
+    println("Σ: \n $Σ")
+    println("abs(_Σ - Σ): \n $(abs.(_Σ - Σ))")
+    println("mean(abs(_Σ - Σ)): \n $(mean(abs.(_Σ - Σ)))")
+
+    @test Σ ≈ _Σ rtol = 0.1
+end
+
+@testset "Two semi-definite constraints" begin
+    # Randomize groundtruth
+    μ1, Σ1 = randn(mat_dim), random_psd_mat(mat_dim)
+    μ2, Σ2 = randn(mat_dim), random_psd_mat(mat_dim)
+
+    # Generate samples
+    samples1 = rand(MvNormal(μ1, Σ1), n_sample)
+    samples2 = rand(MvNormal(μ2, Σ2), n_sample)
+
+    # Objective function
+    function f((x_L1, x_D1, x_L2, x_D2))
+        -loglikelihood(MvNormal(μ1, decompress_symmetric(x_L1, x_D1)), samples1) - 
+            loglikelihood(MvNormal(μ2, decompress_symmetric(x_L2, x_D2)), samples2)
+    end
+    model = Model(f)
+
+    mat_x0 = random_psd_mat(mat_dim)
+    x0 = [mat_x0[Nonconvex.lowertriangind(mat_x0)], diag(mat_x0)]
+    x0 = [copy.(x0)..., copy.(x0)...]
+
+    lbs = [fill(-Inf, length(x0[1])), zeros(length(x0[2]))]
+    lbs = [copy.(lbs)..., copy.(lbs)...]
+
+    ubs = [fill(Inf, length(x0[1])), fill(Inf, length(x0[2]))]
+    ubs = [copy.(ubs)..., copy.(ubs)...]
+
+    addvar!(model, lbs, ubs)
+
+    add_sd_constraint!(model, x -> sd_constraint(x[1:2]))
+    add_sd_constraint!(model, x -> sd_constraint(x[3:4]))
+
+    alg = SDPBarrierAlg(sub_alg=IpoptAlg())
+    options = SDPBarrierOptions(sub_options=IpoptOptions(max_iter=200))
+    result = optimize(model, alg, x0, options = options)
+
+    minimum, minimizer, optimal_ind = result.minimum, result.minimizer, result.optimal_ind
+    _Σ1 = sd_constraint(minimizer[1:2])
+    _Σ2 = sd_constraint(minimizer[3:4])
+
+    println("result: \n $result")
+
+    println("minimum: \n $minimum")
+    println("minimizer: \n $minimizer")
+    println("_Σ1: \n $(_Σ1)")
+    println("_Σ2: \n $(_Σ2)")
+
+    println("Σ1: \n $Σ1")
+    println("Σ2: \n $Σ2")
+    println("abs(_Σ1 - Σ1): \n $(abs.(_Σ1 - Σ1))")
+    println("abs(_Σ2 - Σ2): \n $(abs.(_Σ2 - Σ2))")
+    println("mean(abs(_Σ1 - Σ1)): \n $(mean(abs.(_Σ1 - Σ1)))")
+    println("mean(abs(_Σ2 - Σ2)): \n $(mean(abs.(_Σ2 - Σ2)))")
+
+    @test Σ1 ≈ _Σ1 rtol = 0.1
+    @test Σ2 ≈ _Σ2 rtol = 0.1
+end
+
+@testset "Different algorithm" begin
+    # Randomize groundtruth
+    μ, Σ = randn(mat_dim), random_psd_mat(mat_dim)
+
+    # Generate samples
+    samples = rand(MvNormal(μ, Σ), n_sample)
+
+    # Passing sd_constraint but not using SDPBarrierOptions
+    # Objective function
+    function f((x_L, x_D))
+        -loglikelihood(MvNormal(μ, decompress_symmetric(x_L, x_D)), samples)
+    end
+    model = Model(f)
+
+    mat_x0 = random_psd_mat(mat_dim)
+    x0 = [mat_x0[Nonconvex.lowertriangind(mat_x0)], diag(mat_x0)]
+    lbs = [fill(-Inf, length(x0[1])), zeros(length(x0[2]))]
+    ubs = [fill(Inf, length(x0[1])), fill(Inf, length(x0[2]))]
+    addvar!(model, lbs, ubs)
+    add_sd_constraint!(model, sd_constraint)
+
+    result = optimize(model, IpoptAlg(), x0, options=IpoptOptions(max_iter=1, first_order=true))
+end