In [1]:
## Phase Retrieval example.
# The Phase Retrieval problem seeks to identify
# missing Phase information in imaging problems.
# Classically, it involves a Fourier transform
# and complex variables, but we can study a simple
# variation on the problem.
#
# Given a full-rank matrix A, we are given the
# result of |Ax|.^2 for some hidden vector x,
# (that is, if y = Ax, then we get z_i = |y_i|^2).
# The goal is to the vector x, or a vector that
# has the same norm.

using Random
using LinearAlgebra
Random.seed!(0)
n = 10000
A = randn(n,4)/sqrt(n)
xtrue = [1,0,-1,0]
y = A*xtrue # we don't learn this
z = y.^2 # we only learn this!
# The goal is to get x from z and A!
Out[1]:
10000-element Array{Float64,1}:
 8.731468844193518e-5  
 7.849816728598306e-6  
 1.968606570202451e-5  
 6.653549879680396e-5  
 2.4057369052463896e-5 
 0.00018127177961734864
 4.2158469256264887e-5 
 2.43577848868755e-5   
 2.8034177406233153e-5 
 0.00016317875699489766
 0.00038173091775384535
 0.00029041006054987724
 4.7579876068907034e-5 
 ⋮                     
 2.056941336704841e-5  
 0.0001282448865343001 
 0.0008477646568006551 
 5.294567882390483e-5  
 1.5416847309423263e-5 
 0.00010231563820449386
 0.00012993911483813377
 3.490031021124642e-5  
 1.254712968729905e-6  
 1.913151051310119e-6  
 8.334621053050095e-5  
 0.0008676707498697796
In [2]:
## There is a very good idea that results from using
# the largest eigenvector (this is the eigenvector
# associated with the largest eigenvalue)
# of the matrix A'*Diag(z)*A
#
x = vec(eigen(Symmetric(A'*Diagonal(z)*A)).vectors[:,end])
Out[2]:
4-element Array{Float64,1}:
  0.6991945043822387  
  0.031570773166908435
 -0.7138804075375076  
  0.022474319954127936
In [3]:
##
ntrial = 50
ns = 10.0.^(range(1;length=ntrial,stop=5)) # log-spaced values of n
c = zeros(0)
for n in ns
  n = ceil(Int,n)
  A = randn(n,4)/sqrt(n)
  y = A*xtrue
  z = y.^2
  x = vec(eigen(Symmetric(A'*Diagonal(z)*A)).vectors[:,end])
  @show x
  push!(c, abs(x'*xtrue)/norm(xtrue))
end

syntax: variable "n" declared both local and global

Stacktrace:
 [1] top-level scope at /Users/dgleich/.julia/packages/IJulia/fRegO/src/kernel.jl:52
In [4]:
##
using Plots
pyplot()
myc = c
scatter(ns, myc; xscale=:log10)
#gui()

DivideError: integer division error

Stacktrace:
 [1] div at ./int.jl:228 [inlined]
 [2] fld at ./int.jl:237 [inlined]
 [3] mod at ./int.jl:217 [inlined]
 [4] mod1(::Int64, ::Int64) at ./operators.jl:777
 [5] fix_xy_lengths!(::Plots.Plot{Plots.PyPlotBackend}, ::Plots.Series) at /Users/dgleich/.julia/packages/Plots/h3o4c/src/backends/pyplot.jl:226
 [6] py_add_series(::Plots.Plot{Plots.PyPlotBackend}, ::Plots.Series) at /Users/dgleich/.julia/packages/Plots/h3o4c/src/backends/pyplot.jl:377
 [7] _before_layout_calcs(::Plots.Plot{Plots.PyPlotBackend}) at /Users/dgleich/.julia/packages/Plots/h3o4c/src/backends/pyplot.jl:975
 [8] prepare_output(::Plots.Plot{Plots.PyPlotBackend}) at /Users/dgleich/.julia/packages/Plots/h3o4c/src/plot.jl:254
 [9] show(::Base64.Base64EncodePipe, ::MIME{Symbol("image/png")}, ::Plots.Plot{Plots.PyPlotBackend}) at /Users/dgleich/.julia/packages/Plots/h3o4c/src/output.jl:197
 [10] #base64encode#3(::Nothing, ::typeof(Base64.base64encode), ::Function, ::MIME{Symbol("image/png")}, ::Vararg{Any,N} where N) at /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.2/Base64/src/encode.jl:206
 [11] base64encode(::Function, ::MIME{Symbol("image/png")}, ::Vararg{Any,N} where N) at /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.2/Base64/src/encode.jl:203
 [12] _ijulia_display_dict(::Plots.Plot{Plots.PyPlotBackend}) at /Users/dgleich/.julia/packages/Plots/h3o4c/src/ijulia.jl:50
 [13] display_dict(::Plots.Plot{Plots.PyPlotBackend}) at /Users/dgleich/.julia/packages/Plots/h3o4c/src/init.jl:42
 [14] #invokelatest#1 at ./essentials.jl:790 [inlined]
 [15] invokelatest at ./essentials.jl:789 [inlined]
 [16] execute_request(::ZMQ.Socket, ::IJulia.Msg) at /Users/dgleich/.julia/packages/IJulia/fRegO/src/execute_request.jl:112
 [17] #invokelatest#1 at ./essentials.jl:790 [inlined]
 [18] invokelatest at ./essentials.jl:789 [inlined]
 [19] eventloop(::ZMQ.Socket) at /Users/dgleich/.julia/packages/IJulia/fRegO/src/eventloop.jl:8
 [20] (::getfield(IJulia, Symbol("##15#18")))() at ./task.jl:268