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 039import org.jblas.util.Functions; 040 041/** 042 * Solving linear equations. 043 */ 044public class Solve { 045 /** Solves the linear equation A*X = B. */ 046 public static DoubleMatrix solve(DoubleMatrix A, DoubleMatrix B) { 047 A.assertSquare(); 048 DoubleMatrix X = B.dup(); 049 int[] ipiv = new int[B.rows]; 050 SimpleBlas.gesv(A.dup(), ipiv, X); 051 return X; 052 } 053 054 /** Solves the linear equation A*X = B for symmetric A. */ 055 public static DoubleMatrix solveSymmetric(DoubleMatrix A, DoubleMatrix B) { 056 A.assertSquare(); 057 DoubleMatrix X = B.dup(); 058 int[] ipiv = new int[B.rows]; 059 SimpleBlas.sysv('U', A.dup(), ipiv, X); 060 return X; 061 } 062 063 064 /** Solves the linear equation A*X = B for symmetric and positive definite A. */ 065 public static DoubleMatrix solvePositive(DoubleMatrix A, DoubleMatrix B) { 066 A.assertSquare(); 067 DoubleMatrix X = B.dup(); 068 SimpleBlas.posv('U', A.dup(), X); 069 return X; 070 } 071 072 /** Computes the Least Squares solution for over or underdetermined 073 * linear equations A*X = B 074 * 075 * In the overdetermined case, when m > n, that is, there are more equations than 076 * variables, it computes the least squares solution of X -> ||A*X - B ||_2. 077 * 078 * In the underdetermined case, when m < n (less equations than variables), there are infinitely 079 * many solutions and it computes the minimum norm solution. 080 * 081 * @param A an (m,n) matrix 082 * @param B a (m,k) matrix 083 * @return either the minimum norm or least squares solution. 084 */ 085 public static DoubleMatrix solveLeastSquares(DoubleMatrix A, DoubleMatrix B) { 086 if (B.rows < A.columns) { 087 DoubleMatrix X = DoubleMatrix.concatVertically(B, new DoubleMatrix(A.columns - B.rows, B.columns)); 088 SimpleBlas.gelsd(A.dup(), X); 089 return X; 090 } else { 091 DoubleMatrix X = B.dup(); 092 SimpleBlas.gelsd(A.dup(), X); 093 return X.getRange(0, A.columns, 0, B.columns); 094 } 095 } 096 097 /** 098 * Computes the pseudo-inverse. 099 * 100 * Note, this function uses the solveLeastSquares and might produce different numerical 101 * solutions for the underdetermined case than matlab. 102 * 103 * @param A rectangular matrix 104 * @return matrix P such that A*P*A = A and P*A*P = P. 105 */ 106 public static DoubleMatrix pinv(DoubleMatrix A) { 107 return solveLeastSquares(A, DoubleMatrix.eye(A.rows)); 108 } 109 110//BEGIN 111 // The code below has been automatically generated. 112 // DO NOT EDIT! 113 /** Solves the linear equation A*X = B. */ 114 public static FloatMatrix solve(FloatMatrix A, FloatMatrix B) { 115 A.assertSquare(); 116 FloatMatrix X = B.dup(); 117 int[] ipiv = new int[B.rows]; 118 SimpleBlas.gesv(A.dup(), ipiv, X); 119 return X; 120 } 121 122 /** Solves the linear equation A*X = B for symmetric A. */ 123 public static FloatMatrix solveSymmetric(FloatMatrix A, FloatMatrix B) { 124 A.assertSquare(); 125 FloatMatrix X = B.dup(); 126 int[] ipiv = new int[B.rows]; 127 SimpleBlas.sysv('U', A.dup(), ipiv, X); 128 return X; 129 } 130 131 132 /** Solves the linear equation A*X = B for symmetric and positive definite A. */ 133 public static FloatMatrix solvePositive(FloatMatrix A, FloatMatrix B) { 134 A.assertSquare(); 135 FloatMatrix X = B.dup(); 136 SimpleBlas.posv('U', A.dup(), X); 137 return X; 138 } 139 140 /** Computes the Least Squares solution for over or underdetermined 141 * linear equations A*X = B 142 * 143 * In the overdetermined case, when m > n, that is, there are more equations than 144 * variables, it computes the least squares solution of X -> ||A*X - B ||_2. 145 * 146 * In the underdetermined case, when m < n (less equations than variables), there are infinitely 147 * many solutions and it computes the minimum norm solution. 148 * 149 * @param A an (m,n) matrix 150 * @param B a (m,k) matrix 151 * @return either the minimum norm or least squares solution. 152 */ 153 public static FloatMatrix solveLeastSquares(FloatMatrix A, FloatMatrix B) { 154 if (B.rows < A.columns) { 155 FloatMatrix X = FloatMatrix.concatVertically(B, new FloatMatrix(A.columns - B.rows, B.columns)); 156 SimpleBlas.gelsd(A.dup(), X); 157 return X; 158 } else { 159 FloatMatrix X = B.dup(); 160 SimpleBlas.gelsd(A.dup(), X); 161 return X.getRange(0, A.columns, 0, B.columns); 162 } 163 } 164 165 /** 166 * Computes the pseudo-inverse. 167 * 168 * Note, this function uses the solveLeastSquares and might produce different numerical 169 * solutions for the underdetermined case than matlab. 170 * 171 * @param A rectangular matrix 172 * @return matrix P such that A*P*A = A and P*A*P = P. 173 */ 174 public static FloatMatrix pinv(FloatMatrix A) { 175 return solveLeastSquares(A, FloatMatrix.eye(A.rows)); 176 } 177 178//END 179}