Efficient full-chip thermal simulation is among the most challenging problems facing the EDA industry today, due to the need for solution of very large systems of equations that require unreasonably long computational times. However, in most cases, temperature is not required to be computed at every point of the IC but only at certain hotspots, in order to assess the circuit's compliance with thermal specifications. This makes the thermal analysis problem amenable to Model Order Reduction techniques. System-theoretic techniques like Balanced Truncation offer very reliable bounds for the approximation error, which can be used to control the order and accuracy of the reduced models during creation, at the expense of greater computational complexity to create them. In this paper, we propose a computationally efficient low-rank Balanced Truncation algorithm based on extended Krylov subspace method, which retains all the system-theoretic advantages in the reduction of model order for fast hotspot thermal simulation. Experimental results demonstrate around 97% order reduction and very tight accuracy bounds.