Linear optical circuits of growing complexity are playing an increasing role in emerging photonic quantum technologies. Individual photonic devices are typically described by a unitary matrix containing amplitude and phase information, the characterisation of which is a key task. We present a constructive scheme to retrieve the unitary matrix describing an arbitrary linear optical device using data obtained from one-photon and two-photon ensembles. The scheme is stable on the arbitrarily increasable length scale of the photon packet and independent of photon loss at input and output ports of the device. We find a one-to-one correspondence between ideal data and unitary matrix, and identify the class of non-unitary matrices capable of reproducing the data. The method is extended for coherent state probes, which can simulate two-photon statistics with a reduced visibility. We analyse the performance of reconstruction to simulated noise.