How to return a Symmetric matrix?

General Usage
1

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[1])^2 + b*(x[2] - x[1]^2)^2
end

function h!(H, x; a=1, b=100)
    H[1, 1] = 2 + 8b*x[1]^2 - 4b*(x[2] - x[1]^2)
    H[1, 2] = -4b*x[1]
    #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?

2

return Symmetric(copy(H))?

3

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
4

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

5

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

6

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

7

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

8

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?

9

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[1]^2 - 4b*(x[2] - x[1]^2)
           H.data[1, 2] = -4b*x[1]
           #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
10

are those .data fields present because H is Symmetric?

11

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