# How to return a Symmetric matrix?

I’m constructing the Hessian of the Rosenbrock function. I want the return of `h!` to be of type `Symmetric`.

``````using LinearAlgebra

# for completeness
function f(x; a=1, b=100)
(a - x)^2 + b*(x - x^2)^2
end

function h!(H, x; a=1, b=100)
H[1, 1] = 2 + 8b*x^2 - 4b*(x - x^2)
H[1, 2] = -4b*x
#H[2, 1] = H[1, 2] I want to save this step
H[2, 2] = 2b
return H = copy(Symmetric(H))
end
``````

Now let’s use those functions:

``````julia> H=Matrix{Real}(undef, 2,2)
2×2 Matrix{Real}:
#undef  #undef
#undef  #undef

julia> h!(H, [1.0, 1.0])
2×2 Symmetric{Real, Matrix{Real}}:
802.0  -400.0
-400.0   200

# however
julia> H
2×2 Matrix{Real}:
802.0   -400.0
#undef   200
``````

Changing the last line of `h!` to `copy!(H, Symmetric(H))` returns

``````julia> H
2×2 Matrix{Real}:
802.0  -400.0
-400.0   200
``````

How can I make `h!` return a matrix of type `Symmetric`?

`return Symmetric(copy(H))`?

Your function is returning a `Symmetric` matrix, it’s just not changing the type of `H`. As far as I know you can’t mutate something to a different type. Once you assign a type to something, it’s always going to be that type.

1 Like

Sorry, one more step: you have to assign the result to `H`.

Returns `Symmetric`. It doesn’t change the type of `H`.

But this: `H = Symmetric(H)` changes the type of `H`.

@amrods That rebinds the variable name `H` to a new object but it doesn’t mutate the type of `H` itself.

I see, so I should have started with a symmetric matrix to begin with? For example restricting the type of `H` in as input to `g!` to be `Symmetric`?

Yes I think so, and then you should modify the upper part of the “parent”:

``````julia> H0 = Symmetric(rand(2,2))
2×2 Symmetric{Float64, Matrix{Float64}}:
0.195095  0.336915
0.336915  0.599351

julia> function h!(H::Symmetric, x; a=1, b=100)
H.data[1, 1] = 2 + 8b*x^2 - 4b*(x - x^2)
H.data[1, 2] = -4b*x
#H[2, 1] = H[1, 2] I want to save this step
H.data[2, 2] = 2b
H
end
h! (generic function with 1 method)

julia> h!(H0, ones(2))
2×2 Symmetric{Float64, Matrix{Float64}}:
802.0  -400.0
-400.0   200.0
``````
2 Likes

are those `.data` fields present because `H` is `Symmetric`?

Yes, `S = Symmetric(A)` is a view of the parent `A`. And `A` can be reached with `S.data`. For completeness, `S` also contains a character, either `'u'` or `'l'` to tell if the view uses the upper or lower part of `A`. At the REPL, a quick way to check what’s in an object is to place a dot and hit TAB to autocomplete:

``````julia> H0 = Symmetric(rand(2,2))
2×2 Symmetric{Float64, Matrix{Float64}}:
0.422018  0.115722
0.115722  0.732627

julia> H0. # press TAB twice here with the cursor after the dot
data uplo

julia> H0.data # the non-symmetric parent
2×2 Matrix{Float64}:
0.422018   0.115722
0.0112565  0.732627

julia> H0.uplo # tells that we use the upper part of the parent, i.e., H0[1,2] == H0.parent[1,2]
'U': ASCII/Unicode U+0055 (category Lu: Letter, uppercase)
``````
2 Likes