Large antenna arrays can steer narrow beams towards a target area, and thus improve the communications capacity of wireless channels and the fidelity of radio sensing. Hardware that is capable of continuously-variable phase shifts is expensive, presenting scaling challenges. PIN diodes that apply only discrete phase shifts are promising and cost-effective; however, unlike continuous phase shifters, finding the best phase configuration across elements is an NP-hard optimization problem. Thus, the complexity of optimization becomes a new bottleneck for large-antenna arrays. To address this challenge, this paper suggests a procedure for converting the optimization objective function from a ratio of quadratic functions to a sequence of more easily solvable quadratic unconstrained binary optimization (QUBO) sub-problems. This conversion is an exact equivalence, and the resulting QUBO forms are standard input formats for various physics-inspired optimization methods. We demonstrate that a simulated annealing approach is very effective for solving these sub-problems, and we give performance metrics for several large array types optimized by this technique. Through numerical experiments, we report 3D beamforming performance for extra-large arrays with up to 10,000 elements.