001// --- BEGIN LICENSE BLOCK --- 002/* 003 * Copyright (c) 2009, Mikio L. Braun 004 * All rights reserved. 005 * 006 * Redistribution and use in source and binary forms, with or without 007 * modification, are permitted provided that the following conditions are 008 * met: 009 * 010 * * Redistributions of source code must retain the above copyright 011 * notice, this list of conditions and the following disclaimer. 012 * 013 * * Redistributions in binary form must reproduce the above 014 * copyright notice, this list of conditions and the following 015 * disclaimer in the documentation and/or other materials provided 016 * with the distribution. 017 * 018 * * Neither the name of the Technische Universität Berlin nor the 019 * names of its contributors may be used to endorse or promote 020 * products derived from this software without specific prior 021 * written permission. 022 * 023 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 024 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 025 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 026 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 027 * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 028 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 029 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 030 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 031 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 032 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 033 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 034 */ 035// --- END LICENSE BLOCK --- 036 037package org.jblas; 038 039/** 040 * Solving linear equations. 041 */ 042public class Solve { 043 /** Solves the linear equation A*X = B. */ 044 public static DoubleMatrix solve(DoubleMatrix A, DoubleMatrix B) { 045 A.assertSquare(); 046 DoubleMatrix X = B.dup(); 047 int[] ipiv = new int[B.rows]; 048 SimpleBlas.gesv(A.dup(), ipiv, X); 049 return X; 050 } 051 052 /** Solves the linear equation A*X = B for symmetric A. */ 053 public static DoubleMatrix solveSymmetric(DoubleMatrix A, DoubleMatrix B) { 054 A.assertSquare(); 055 DoubleMatrix X = B.dup(); 056 int[] ipiv = new int[B.rows]; 057 SimpleBlas.sysv('U', A.dup(), ipiv, X); 058 return X; 059 } 060 061 062 /** Solves the linear equation A*X = B for symmetric and positive definite A. */ 063 public static DoubleMatrix solvePositive(DoubleMatrix A, DoubleMatrix B) { 064 A.assertSquare(); 065 DoubleMatrix X = B.dup(); 066 SimpleBlas.posv('U', A.dup(), X); 067 return X; 068 } 069 070 /** Computes the Least Squares solution for over or underdetermined 071 * linear equations A*X = B 072 * 073 * In the overdetermined case, when m > n, that is, there are more equations than 074 * variables, it computes the least squares solution of X -> ||A*X - B ||_2. 075 * 076 * In the underdetermined case, when m < n (less equations than variables), there are infinitely 077 * many solutions and it computes the minimum norm solution. 078 * 079 * @param A an (m,n) matrix 080 * @param B a (m,k) matrix 081 * @return either the minimum norm or least squares solution. 082 */ 083 public static DoubleMatrix solveLeastSquares(DoubleMatrix A, DoubleMatrix B) { 084 if (B.rows < A.columns) { 085 DoubleMatrix X = DoubleMatrix.concatVertically(B, new DoubleMatrix(A.columns - B.rows, B.columns)); 086 SimpleBlas.gelsd(A.dup(), X); 087 return X; 088 } else { 089 DoubleMatrix X = B.dup(); 090 SimpleBlas.gelsd(A.dup(), X); 091 return X.getRange(0, A.columns, 0, B.columns); 092 } 093 } 094 095 /** 096 * Computes the pseudo-inverse. 097 * 098 * Note, this function uses the solveLeastSquares and might produce different numerical 099 * solutions for the underdetermined case than matlab. 100 * 101 * @param A rectangular matrix 102 * @return matrix P such that A*P*A = A and P*A*P = P. 103 */ 104 public static DoubleMatrix pinv(DoubleMatrix A) { 105 return solveLeastSquares(A, DoubleMatrix.eye(A.rows)); 106 } 107 108//BEGIN 109 // The code below has been automatically generated. 110 // DO NOT EDIT! 111 /** Solves the linear equation A*X = B. */ 112 public static FloatMatrix solve(FloatMatrix A, FloatMatrix B) { 113 A.assertSquare(); 114 FloatMatrix X = B.dup(); 115 int[] ipiv = new int[B.rows]; 116 SimpleBlas.gesv(A.dup(), ipiv, X); 117 return X; 118 } 119 120 /** Solves the linear equation A*X = B for symmetric A. */ 121 public static FloatMatrix solveSymmetric(FloatMatrix A, FloatMatrix B) { 122 A.assertSquare(); 123 FloatMatrix X = B.dup(); 124 int[] ipiv = new int[B.rows]; 125 SimpleBlas.sysv('U', A.dup(), ipiv, X); 126 return X; 127 } 128 129 130 /** Solves the linear equation A*X = B for symmetric and positive definite A. */ 131 public static FloatMatrix solvePositive(FloatMatrix A, FloatMatrix B) { 132 A.assertSquare(); 133 FloatMatrix X = B.dup(); 134 SimpleBlas.posv('U', A.dup(), X); 135 return X; 136 } 137 138 /** Computes the Least Squares solution for over or underdetermined 139 * linear equations A*X = B 140 * 141 * In the overdetermined case, when m > n, that is, there are more equations than 142 * variables, it computes the least squares solution of X -> ||A*X - B ||_2. 143 * 144 * In the underdetermined case, when m < n (less equations than variables), there are infinitely 145 * many solutions and it computes the minimum norm solution. 146 * 147 * @param A an (m,n) matrix 148 * @param B a (m,k) matrix 149 * @return either the minimum norm or least squares solution. 150 */ 151 public static FloatMatrix solveLeastSquares(FloatMatrix A, FloatMatrix B) { 152 if (B.rows < A.columns) { 153 FloatMatrix X = FloatMatrix.concatVertically(B, new FloatMatrix(A.columns - B.rows, B.columns)); 154 SimpleBlas.gelsd(A.dup(), X); 155 return X; 156 } else { 157 FloatMatrix X = B.dup(); 158 SimpleBlas.gelsd(A.dup(), X); 159 return X.getRange(0, A.columns, 0, B.columns); 160 } 161 } 162 163 /** 164 * Computes the pseudo-inverse. 165 * 166 * Note, this function uses the solveLeastSquares and might produce different numerical 167 * solutions for the underdetermined case than matlab. 168 * 169 * @param A rectangular matrix 170 * @return matrix P such that A*P*A = A and P*A*P = P. 171 */ 172 public static FloatMatrix pinv(FloatMatrix A) { 173 return solveLeastSquares(A, FloatMatrix.eye(A.rows)); 174 } 175 176//END 177}