Hamiltonian diagonalization is at the heart of understanding physical properties and practical applications of quantum systems. It is highly desired to design quantum algorithms that can speedup Hamiltonian diagonalization, especially those can be implemented on near-term quantum devices. In this work, we propose a variational quantum algorithm for Hamiltonian diagonalization (VQHD) of quantum systems, which explores the important physical properties, such as temperature, locality, and correlation, of the system. The key idea is that the thermal states of the system encode the information of eigenvalues and eigenstates of the system Hamiltonian. To obtain the full spectrum of the Hamiltonian, we use a quantum imaginary time evolution algorithm with high temperature, which prepares a thermal state with a small correlation length. With Trotterization, this then allows us to implement each step of imaginary time evolution by a local unitary transformation on only a small number of sites. Diagonalizing these thermal states hence leads to a full knowledge of the Hamiltonian eigensystem. We apply our algorithm to diagonalize local Hamiltonians and return results with high precision. Our VQHD algorithm sheds new light on the applications of near-term quantum computers.