find_K.jl

function build_A(κ) 
    N = length(κ)
    A = zeros(N,N)
    for n=2:N-2
        # A[n,n] = 2*κ[n]
        # A[n,n+1] = 1/4*(κ[n-1] - κ[n+1] - 4κ[n])
        # A[n+1,n] = 1/4*(κ[n+2] - κ[n] - 4κ[n+1])

        #Reference case: we know that AK is symmetric with K = diag(1/κ)
        A[n,n] = -2*(κ[n+1] + κ[n-1])
        A[n,n+1] = κ[n+1]
        A[n+1,n] = κ[n]

    end
    A
end

#attempts to find a K diagonal such that AK is symmetric
#this assumes A has distinct eigenvalues
function find_K(A)
    # For A = V D V^-1, we have A = V D V^T (V V^T)^-1
    # More generally this is true replacing V by V S with S diagonal (changing the scaling of the eigenvectors)
    # Therefore take V S^2 V', where S diagonal is chosen to make V S^2 V' diagonal

    # S satisfies the overdetermined system \sum_k S_k^2 V_ik V_jk = 0 ∀ i ≠ j, so we form that matrix and compute S as its nullspace
    N = size(A,1)
    V = eig(A)[2]
    mat = zeros(div(N*(N-1),2), N)
    for k=1:N
        mat[:,k] = [V[i,k]*V[j,k] for i=1:N, j=1:N if i < j]
    end
    # S^2 should satisfy mat*S2 = 0
    S2 = nullspace(mat)
    @assert size(S2,2) == 1 #hope for unique solution
    S2 = S2[:,1]
    # flip the sign to make it positive if needed
    if S2[1] < 0
        S2 = -S2
    end
    # if those println are 1e-15 then the procedure has succeeded
    println(norm(mat*S2))
    K = V*diagm(S2[:,1])*V'
    Asym = A*K-K*A'
    println(norm(Asym))
    println(norm(K - diagm(diag(K))))
    diag(K)
end
N = 20
κ = cos.(0.1*(1:N)/N)
κ = 1+rand(N)
# κ = ones(n)
κdiag = diagm(κ)
A = build_A(κ)

trim = 4
A = A[trim:N-trim,trim:N-trim]
κ = κ[trim:N-trim]

K = find_K(A)
κ.*K #κ.*K should be constant (K = 1/κ up to scaling) for the reference case